Hydrological Analysis Report

1. DATA PREPARATION

In this section, we import the hydrological data and perform initial processing, which includes checking for any missing values and cleaning the dataset. If missing data is found, appropriate methods will be used to address them.

# Import dataset from excel file
main_data = read_excel("Savinja-VelikoSirje.xlsx")

# Create dataframe
main_data = as.data.frame(main_data)

# Convert the 'Date' column to Date type
main_data$Date <- as.Date(main_data$Date, format = "%d.%m.%Y")

head(main_data)

##         Date Year Month Day Waterlevel Discharge Watertemp P1_Luce P2_Solcava
## 1 2010-01-01 2010     1   1        251    59.058       6.9     7.0        4.8
## 2 2010-01-02 2010     1   2        246    53.925       7.1     3.0        0.6
## 3 2010-01-03 2010     1   3        239    47.681       6.6     0.1        2.7
## 4 2010-01-04 2010     1   4        231    41.040       5.3     0.0        0.0
## 5 2010-01-05 2010     1   5        228    38.694       4.4     3.3        2.0
## 6 2010-01-06 2010     1   6        226    36.966       4.2     6.2        4.6
##   P3_Lasko P4_GornjiGrad P5_Radegunda Airtemp Evapotrans
## 1      4.3           4.6          2.2     5.7        0.3
## 2      2.1           3.2          2.0     3.5        0.9
## 3      1.5           0.0          0.1    -2.0        0.5
## 4      0.0           0.0          0.0    -3.0        0.2
## 5      7.0           3.1          4.0    -1.8        0.2
## 6      7.0           9.2          8.2    -1.0        0.2

# Structure of dataframe
str(main_data)

## 'data.frame':    4748 obs. of  14 variables:
##  $ Date         : Date, format: "2010-01-01" "2010-01-02" ...
##  $ Year         : num  2010 2010 2010 2010 2010 2010 2010 2010 2010 2010 ...
##  $ Month        : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ Day          : num  1 2 3 4 5 6 7 8 9 10 ...
##  $ Waterlevel   : num  251 246 239 231 228 226 223 220 223 223 ...
##  $ Discharge    : num  59.1 53.9 47.7 41 38.7 ...
##  $ Watertemp    : num  6.9 7.1 6.6 5.3 4.4 4.2 4.2 4.2 3.7 3.6 ...
##  $ P1_Luce      : num  7 3 0.1 0 3.3 6.2 2.5 2 37.6 2.4 ...
##  $ P2_Solcava   : num  4.8 0.6 2.7 0 2 4.6 0.4 2 32.2 4.5 ...
##  $ P3_Lasko     : num  4.3 2.1 1.5 0 7 7 8.1 3.1 20.1 3.2 ...
##  $ P4_GornjiGrad: num  4.6 3.2 0 0 3.1 9.2 3.9 2.8 36.1 4.3 ...
##  $ P5_Radegunda : num  2.2 2 0.1 0 4 8.2 2.8 1.5 34.1 3.8 ...
##  $ Airtemp      : num  5.7 3.5 -2 -3 -1.8 -1 -0.3 -0.4 0.6 0.9 ...
##  $ Evapotrans   : num  0.3 0.9 0.5 0.2 0.2 0.2 0.2 0.2 0.2 0.2 ...

# Check for missing values
missing_data <- colSums(is.na(main_data))
missing_data

##          Date          Year         Month           Day    Waterlevel 
##             0             0             0             0           183 
##     Discharge     Watertemp       P1_Luce    P2_Solcava      P3_Lasko 
##             0           224             0             0           184 
## P4_GornjiGrad  P5_Radegunda       Airtemp    Evapotrans 
##             0             0             0             0

1.1. Fill Missing Waterlevel Values

sum(is.na(main_data$Waterlevel))

## [1] 183

# Fit linear regression model using Discharge
model <- lm(Waterlevel ~ Discharge, data = main_data, na.action = na.exclude)
model

## 
## Call:
## lm(formula = Waterlevel ~ Discharge, data = main_data, na.action = na.exclude)
## 
## Coefficients:
## (Intercept)    Discharge  
##      173.35         0.88

# Predict missing Waterlevel
missing_waterlevel <- which(is.na(main_data$Waterlevel))
predicted_values <- predict(model, newdata = main_data[missing_waterlevel, ])
predicted_values

##     2009     2010     2011     2012     2013     2014     2015     2016 
## 197.4418 195.2383 193.6437 192.6052 191.6953 190.5953 189.9546 217.9093 
##     2017     2018     2019     2020     2021     2022     2023     2024 
## 203.1065 195.1661 192.8129 191.1277 190.0981 188.8528 188.0775 187.6472 
##     2025     2026     2027     2028     2029     2030     2031     2032 
## 186.7047 186.0421 185.6117 185.0608 185.3196 184.5997 188.0907 201.2136 
##     2033     2034     2035     2036     2037     2038     2039     2040 
## 197.9663 196.5495 200.5580 194.4436 227.3713 227.5033 203.9716 197.2421 
##     2041     2042     2043     2044     2045     2046     2047     2048 
## 196.7756 195.4697 192.1256 190.0514 188.7631 187.7326 186.7461 186.1406 
##     2049     2050     2051     2052     2053     2054     2055     2056 
## 185.2861 184.8743 184.1685 183.8543 183.8367 184.4607 186.8376 186.7883 
##     2057     2058     2059     2060     2061     2062     2063     2064 
## 187.9411 205.1640 187.8223 185.2923 184.1333 183.6141 183.5384 187.5610 
##     2065     2066     2067     2068     2069     2070     2071     2072 
## 186.4134 184.5117 183.6141 183.0482 182.6170 182.1700 181.9465 182.0213 
##     2073     2074     2075     2076     2077     2078     2079     2080 
## 184.0268 222.2223 237.3929 208.2124 196.2529 191.6548 189.4240 187.9271 
##     2081     2082     2083     2084     2085     2086     2087     2088 
## 186.8904 186.0421 185.4190 198.7214 193.9270 190.3366 188.5862 187.6692 
##     2089     2090     2091     2092     2093     2094     2095     2096 
## 186.9133 185.9030 185.3433 184.9491 193.3603 218.6925 220.0477 202.5609 
##     2097     2098     2099     2100     2101     2102     2103     2104 
## 195.5586 191.9012 190.0285 188.8308 187.8611 187.3700 186.7866 187.7150 
##     2105     2106     2107     2108     2109     2110     2111     2112 
## 186.8904 229.9902 225.4678 207.6984 198.9590 264.2088 261.3857 228.6218 
##     2113     2114     2115     2116     2117     2118     2119     2120 
## 307.3579 510.5842 405.4642 279.2572 245.3932 254.7803 258.4201 232.4146 
##     2121     2122     2123     2124     2125     2126     2127     2128 
## 220.7342 213.7688 208.5626 205.1464 202.5169 200.2508 198.5858 197.0643 
##     2129     2130     2131     2132     2133     2134     2135     2136 
## 195.4847 194.1233 192.7856 191.8203 191.2755 190.7167 190.0972 189.4900 
##     2137     2138     2139     2140     2141     2142     2143     2144 
## 189.0104 188.4067 188.0775 187.7106 187.2635 186.8490 186.4812 186.2198 
##     2145     2146     2147     2148     2149     2150     2151     2152 
## 185.8045 185.8045 185.8837 185.1559 184.9508 185.6901 218.8606 224.1381 
##     2153     2154     2155     2156     2157     2158     2159     2160 
## 207.3922 198.0130 195.1529 193.7554 191.7763 190.3814 189.4451 188.7842 
##     2161     2162     2163     2164     2165     2166     2167     2168 
## 188.8291 188.8520 189.2603 188.8317 187.5407 186.9106 186.4011 186.2401 
##     2169     2170     2171     2172     2173     2174     2175     2176 
## 185.8247 185.6311 184.9667 184.8558 184.6728 184.2935 184.0268 183.8367 
##     2177     2178     2179     2180     2181     2182     2183     2184 
## 183.6123 183.6123 183.6123 183.6123 183.5454 183.5789 183.2603 183.0808 
##     2185     2186     2187     2188     2189     2190     2191 
## 182.8370 182.8212 182.8212 182.8212 182.8212 182.8212 182.5871

# Fill missing Waterlevel values
main_data$Waterlevel[missing_waterlevel] <- predicted_values

# Check Waterlevel again
sum(is.na(main_data$Waterlevel))

## [1] 0

1.2. Fill Missing Watertemp Values

missing_watertemps = which(is.na(main_data$Watertemp))
missing_watertemps

##   [1] 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992
##  [16] 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007
##  [31] 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022
##  [46] 2023 2024 2025 2026 2027 2028 2029 2030 2031 2032 2033 2034 2035 2036 2037
##  [61] 2038 2039 2040 2041 2042 2043 2044 2045 2046 2047 2048 2049 2050 2051 2052
##  [76] 2053 2054 2055 2056 2057 2058 2059 2060 2061 2062 2063 2064 2065 2066 2067
##  [91] 2068 2069 2070 2071 2072 2073 2074 2075 2076 2077 2078 2079 2080 2081 2082
## [106] 2083 2084 2085 2086 2087 2088 2089 2090 2091 2092 2093 2094 2095 2096 2097
## [121] 2098 2099 2100 2101 2102 2103 2104 2105 2106 2107 2108 2109 2110 2111 2112
## [136] 2113 2114 2115 2116 2117 2118 2119 2120 2121 2122 2123 2124 2125 2126 2127
## [151] 2128 2129 2130 2131 2132 2133 2134 2135 2136 2137 2138 2139 2140 2141 2142
## [166] 2143 2144 2145 2146 2147 2148 2149 2150 2151 2152 2153 2154 2155 2156 2157
## [181] 2158 2159 2160 2161 2162 2163 2164 2165 2166 2167 2168 2169 2170 2171 2172
## [196] 2173 2174 2175 2176 2177 2178 2179 2180 2181 2182 2183 2184 2185 2186 2187
## [211] 2188 2189 2190 2191 2192 2193 2194 2195 2196 2197 2198 2199 2200 2201

# Loop missing watertemp and fill 
for (i in missing_watertemps){
  current_year <- main_data$Year[i]
  date_without_year <- format(main_data$Date[i], "%m-%d")
  
  # Find Watertemp and Airtemp for same date in previous years
  previous_data <- main_data %>% 
    filter(format(Date, "%m-%d") == date_without_year & Year < current_year) %>% 
    dplyr::select(Watertemp, Airtemp)

  
  # Calculate mean of Watertemp and Airtemp
  if(nrow(previous_data) > 0){
    watertemp_mean <- mean(previous_data$Watertemp)
    airtemp_mean <- mean(previous_data$Airtemp)
    
    main_data$Watertemp[i] <- watertemp_mean + (main_data$Airtemp[i] - airtemp_mean)
  }
}

# Check Watertemp
sum(is.na(main_data$Watertemp))

## [1] 0

1.3. Fill Missing Values for P3_Lasko

missing_lasko <-  which(is.na(main_data$P3_Lasko))

# Fill missing values with average of other stations
for (i in missing_lasko){
  other_prec <- main_data[i, c("P1_Luce", "P2_Solcava", "P4_GornjiGrad", "P5_Radegunda")]
  main_data$P3_Lasko[i] <- mean(as.numeric(other_prec))
}

# Check P3_Lasko
sum(is.na(main_data$P3_Lasko))

## [1] 0

# Check main_data for any missing values
colSums(is.na(main_data))

##          Date          Year         Month           Day    Waterlevel 
##             0             0             0             0             0 
##     Discharge     Watertemp       P1_Luce    P2_Solcava      P3_Lasko 
##             0             0             0             0             0 
## P4_GornjiGrad  P5_Radegunda       Airtemp    Evapotrans 
##             0             0             0             0

2. BASIC ANALYSIS

In this section, we perform basic hydrological analysis. This includes calculating descriptive statistics for key variables such as discharge, precipitation, evapotranspiration, and air temperature. We will also plot graphs for these variables over the entire period and for the year 2022. Additionally, we will calculate correlations between precipitation stations, water balance, runoff coefficient, and the ratio between evapotranspiration and precipitation.

2.1. Calculate Descriptive Statistics

# Calculate descriptive statistics for Discharge
discharge_stats <- data.frame(
  Statistic = c("Mean", "Standard Deviation", "Minimum", "Maximum", "Median"),
  Values = c(
    mean(main_data$Discharge),
    sd(main_data$Discharge),
    min(main_data$Discharge),
    max(main_data$Discharge),
    median(main_data$Discharge))
) 

discharge_stats

##            Statistic    Values
## 1               Mean  38.48014
## 2 Standard Deviation  51.68336
## 3            Minimum   5.51500
## 4            Maximum 859.74400
## 5             Median  22.03150

# Summary statistics for Discharge
summary(main_data$Discharge)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   5.515  13.966  22.032  38.480  40.252 859.744

# Calculate descriptive statistics for Precipitation 
precipitation_stats <- data.frame(
  Statistic = c("Mean", "Standard Deviation", "Minimum", "Maximum", "Median"),
  Values = c(
    mean(c(main_data$P1_Luce, main_data$P2_Solcava, main_data$P3_Lasko, main_data$P4_GornjiGrad, main_data$P5_Radegunda)),
    sd(c(main_data$P1_Luce, main_data$P2_Solcava, main_data$P3_Lasko, main_data$P4_GornjiGrad, main_data$P5_Radegunda)),
    min(c(main_data$P1_Luce, main_data$P2_Solcava, main_data$P3_Lasko, main_data$P4_GornjiGrad, main_data$P5_Radegunda)),
    max(c(main_data$P1_Luce, main_data$P2_Solcava, main_data$P3_Lasko, main_data$P4_GornjiGrad, main_data$P5_Radegunda)),
    median(c(main_data$P1_Luce, main_data$P2_Solcava, main_data$P3_Lasko, main_data$P4_GornjiGrad, main_data$P5_Radegunda))
  )
)

precipitation_stats

##            Statistic     Values
## 1               Mean   3.925959
## 2 Standard Deviation   9.844119
## 3            Minimum   0.000000
## 4            Maximum 135.400000
## 5             Median   0.000000

# Summary statistics for Precipitation
summary(c(main_data$P1_Luce, main_data$p2_Solcava, main_data$P3_Lasko, main_data$P4_GornjiGrad, main_data$P5_Radegunda))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    0.00    0.00    3.85    2.10  135.40

# P1_Luce
mean(main_data$P1_Luce)

## [1] 4.235573

sd(main_data$P1_Luce)

## [1] 10.816

min(main_data$P1_Luce)

## [1] 0

max(main_data$P1_Luce)

## [1] 116.1

median(main_data$P1_Luce)

## [1] 0

# P2_Solcava
mean(main_data$P2_Solcava)

## [1] 4.23155

sd(main_data$P2_Solcava)

## [1] 10.60273

min(main_data$P2_Solcava)

## [1] 0

max(main_data$P2_Solcava)

## [1] 109.5

median(main_data$P2_Solcava)

## [1] 0

# P3_Lasko
mean(main_data$P3_Lasko)

## [1] 3.255281

sd(main_data$P3_Lasko)

## [1] 8.140033

min(main_data$P3_Lasko)

## [1] 0

max(main_data$P3_Lasko)

## [1] 106.7

median(main_data$P3_Lasko)

## [1] 0

# P4_GornjiGrad
mean(main_data$P4_GornjiGrad)

## [1] 4.119756

sd(main_data$P4_GornjiGrad)

## [1] 10.10802

min(main_data$P4_GornjiGrad)

## [1] 0

max(main_data$P4_GornjiGrad)

## [1] 135.4

median(main_data$P4_GornjiGrad)

## [1] 0

# P5_Radegunda
mean(main_data$P5_Radegunda)

## [1] 3.787637

sd(main_data$P5_Radegunda)

## [1] 9.278106

min(main_data$P5_Radegunda)

## [1] 0

max(main_data$P5_Radegunda)

## [1] 110.5

median(main_data$P5_Radegunda)

## [1] 0

# Calculate descriptive statistics for Air Temperature
airtemp_stats <- data.frame(
  Statistic = c("Mean", "Standard Deviation", "Minimum", "Maximum", "Median"),
  Values = c(
    mean(main_data$Airtemp),
    sd(main_data$Airtemp),
    min(main_data$Airtemp),
    max(main_data$Airtemp),
    median(main_data$Airtemp)
  )
)

airtemp_stats

##            Statistic     Values
## 1               Mean  10.878265
## 2 Standard Deviation   8.140201
## 3            Minimum -13.800000
## 4            Maximum  28.900000
## 5             Median  11.300000

# Summary statistics for Air Temperature
summary(main_data$Airtemp)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -13.80    4.40   11.30   10.88   17.70   28.90

# Calculate descriptive statistics for Evapotranspiration
evapotrans_stats <- data.frame(
  Statistic = c("Mean", "Standard Deviation", "Minimum", "Maximum", "Median"),
  Values = c(
    mean(main_data$Evapotrans),
    sd(main_data$Evapotrans),
    min(main_data$Evapotrans),
    max(main_data$Evapotrans),
    median(main_data$Evapotrans)
  )
)

evapotrans_stats

##            Statistic   Values
## 1               Mean 2.264006
## 2 Standard Deviation 1.713055
## 3            Minimum 0.000000
## 4            Maximum 7.700000
## 5             Median 1.800000

# Summary statistics for Evapotranspiration
summary(main_data$Evapotrans)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.700   1.800   2.264   3.600   7.700

2.2. Plot Graphs for Key Variables

# Plot Discharge over the whole period
plot(main_data$Date, main_data$Discharge, type = "l",
     xlab = "Date", ylab = "Discharge (m^3/s)", main = "Discharge Over Whole Period")

# Plot Discharge for year 2022
data_2022 <- main_data[format(main_data$Date, "%Y") == "2022",]
plot(data_2022$Date, data_2022$Discharge, type = "l",
     xlab = "Date", ylab = "Discharge (m^3/s)", main = "Discharge For Year 2022")

# Plot Precipitation over whole period for all stations
plot(main_data$Date, main_data$P1_Luce, type = "l", col = "blue",
     xlab = "Date", ylab = "Precipitation (mm)", main = "Precipitation Over Whole Period")

# All lines for other stations
lines(main_data$Date, main_data$P2_Solcava, col = "red")
lines(main_data$Date, main_data$P3_Lasko, col = "green")
lines(main_data$Date, main_data$P4_GornjiGrad, col = "purple")
lines(main_data$Date, main_data$P5_Radegunda, col = "orange")

# Add legend
legend("topright", legend = c("P1_Luce", "P2_Solcava", "P3_Lasko", "P4_GornjiGrad", "P5_Radegunda"),
       col = c("blue", "red", "green", "purple", "orange"), lty = 1)

# Plot Precipitation for year 2022 
plot(data_2022$Date, data_2022$P1_Luce, type = "l", col = "blue",
     xlab = "Date", ylab = "Precipitation (mm)", main = "Precipitation for Year 2022")

# Add lines for other stations
lines(data_2022$Date, data_2022$p2_Solcava, col = "red")
lines(data_2022$Date, data_2022$P3_Lasko, col = "green")
lines(data_2022$Date, data_2022$P4_GornjiGrad, col = "purple")
lines(data_2022$Date, data_2022$P5_Radegunda, col = "orange")

# Add a legend
legend("topright", legend = c("P1_Luce", "P2_Solcava", "P3_Lasko", "P4_GornjiGrad", "P5_Radegunda"),
       col = c("blue", "red", "green", "purple", "orange"), lty = 1)

# Calculate the average precipitation across all stations
main_data$Avg_Precipitation <- rowMeans(main_data[, c("P1_Luce", "P2_Solcava", "P3_Lasko", "P4_GornjiGrad", "P5_Radegunda")])

# Plot Average Precipitation over the whole period
plot(main_data$Date, main_data$Avg_Precipitation, type = "l",
     xlab = "Date", ylab = "Average Precipitation (mm)", main = "Average Precipitation Over Whole Period")

# Plot Average Precipitation for year 2022
data_2022 <- main_data[format(main_data$Date, "%Y") == "2022",]
plot(data_2022$Date, data_2022$Avg_Precipitation, type = "l",
     xlab = "Date", ylab = "Average Precipitation (mm)", main = "Average Precipitation For Year 2022")

# Plot Evapotranspiration over whole period
plot(main_data$Date, main_data$Evapotrans, type = "l", col = "blue",
     xlab = "Date", ylab = "Evapotranspiration (mm)", main = "Evapotranspiration Over Whole Period")

# Plot Evapotranspiration for year 2022
data_2022 <- main_data[format(main_data$Date, "%Y") == "2022",]
plot(data_2022$Date, data_2022$Evapotrans, type = "l", col = "blue",
     xlab = "Date", ylab = "Evapotranspiration (mm)", main = "Evapotranspiration For Year 2022")

# Plot Air Temperature over the whole period
plot(main_data$Date, main_data$Airtemp, type = "l", col = "red",
     xlab = "Date", ylab = "Air Temperature (°C)", main = "Air Temperature Over Whole Period")

# Plot Air Temperature for year 2022
plot(data_2022$Date, data_2022$Airtemp, type = "l", col = "red",
     xlab = "Date", ylab = "Air Temperature (°C)", main = "Air Temperature For Year 2022")

2.3. Calculate Correlation Between Precipitation Stations

# Precipitation columns
precipitation_data <- main_data[, c("P1_Luce", "P2_Solcava", "P3_Lasko", "P4_GornjiGrad", "P5_Radegunda")]

# Calculate correlation matrix
cor_matrix <- cor(precipitation_data)

cor_matrix

##                 P1_Luce P2_Solcava  P3_Lasko P4_GornjiGrad P5_Radegunda
## P1_Luce       1.0000000  0.8974698 0.7280284     0.9175150    0.8899014
## P2_Solcava    0.8974698  1.0000000 0.7226416     0.8618554    0.8479587
## P3_Lasko      0.7280284  0.7226416 1.0000000     0.7861900    0.8029485
## P4_GornjiGrad 0.9175150  0.8618554 0.7861900     1.0000000    0.9092337
## P5_Radegunda  0.8899014  0.8479587 0.8029485     0.9092337    1.0000000

# Plot correlogram
corrplot(cor_matrix, method = "shade", type = "full", 
         tl.col = "orange", tl.srt = 45,
         title = "Correlogram of Precipitation Stations", 
          mar = c(0, 0, 2, 0)
)

The correlation analysis for selected precipitation stations in the Savinja River basin reveals high inter-station correlations, with values ranging from 0.72 to 0.92. The strongest correlation is observed between P1_Luce and P4_GornjiGrad (0.92), indicating similar precipitation patterns, likely influenced by geographic or climatic factors. The weakest correlation (0.72) is between P2_Solcava and P3_Lasko, which may reflect localized variations in precipitation. The correlogram visually highlights these relationships, with darker shades representing stronger correlations. Stations such as P1_Luce, P4_GornjiGrad, and P5_Radegunda exhibit particularly high inter-correlations, suggesting potential redundancy in data collection across these locations.

2.4. Calculate Water Balance

# Catchment area in m²
catchment_area <- 1847000000

# Calculate Runoff
main_data$Runoff <- (main_data$Discharge * 3600 * 24 * 1000) / catchment_area

main_data$Discharge_mm <- main_data$Discharge * 3600 * 24 * 1000

# Calculate the water balance
main_data$Water_Balance <- main_data$Avg_Precipitation - main_data$Evapotrans - main_data$Runoff

# Calculated Runoff
main_data$Calculated_Runoff <- main_data$Avg_Precipitation - main_data$Evapotrans

# Plot measured and calculated runoff
yearly_runoff <- main_data %>% 
  group_by(Year) %>% 
  summarise(
    Total_Measured = sum(Runoff),
    Total_Calculated = sum(Calculated_Runoff))

plot(yearly_runoff$Year, yearly_runoff$Total_Measured, 
     xlab = "Date", ylab ="Runoff", type = "l", col = "red", main = "Measured and Calculated Runoff")
lines(yearly_runoff$Year, yearly_runoff$Total_Calculated, lty = 2, col = "green")
legend("topright", legend = c("Measured", "Calculated"), col = c("red", "green"), lty = c(1, 2))

The comparison between measured and calculated runoff for the Savinja River basin shows a general agreement in the annual trends. While the measured runoff (red line) closely aligns with the calculated runoff (green dashed line), slight deviations are observed in certain years, such as 2012 and 2020, which may be attributed to discrepancies in precipitation or evapotranspiration data. Overall, the calculated runoff provides a reasonable approximation of the measured values, validating the water balance approach.

2.5. Calculate Runoff Coefficient

# Group data by Year and calculate total runoff and total precipitation
annual_runoff <- main_data %>%
  group_by(Year) %>%
  summarise(
    Total_Runoff = sum(Runoff),             
    Total_Precipitation = sum(Avg_Precipitation)
  )

# Calculate the runoff coefficient for each year
annual_runoff$Runoff_Coefficient <- annual_runoff$Total_Runoff / annual_runoff$Total_Precipitation

annual_runoff

## # A tibble: 13 × 4
##     Year Total_Runoff Total_Precipitation Runoff_Coefficient
##    <dbl>        <dbl>               <dbl>              <dbl>
##  1  2010         816.               1603.              0.509
##  2  2011         369.               1168.              0.316
##  3  2012         625.               1514.              0.413
##  4  2013         793.               1484.              0.535
##  5  2014        1037.               1892.              0.548
##  6  2015         546.               1205.              0.454
##  7  2016         685.               1435.              0.478
##  8  2017         632.               1514.              0.417
##  9  2018         687.               1314.              0.523
## 10  2019         665.               1512.              0.440
## 11  2020         551.               1445.              0.381
## 12  2021         664.               1403.              0.474
## 13  2022         476.               1151.              0.413

# Plot the runoff coefficient over the years
plot(as.numeric(annual_runoff$Year), annual_runoff$Runoff_Coefficient, type = "b", col = "blue",
     xlab = "Year", ylab = "Runoff Coefficient", main = "Runoff Coefficient Values")

The runoff coefficient values over the years show fluctuations between 0.31 and 0.55. Higher runoff coefficients, such as in 2013 and 2014, indicate a larger proportion of precipitation contributing to runoff, potentially due to higher precipitation intensity or reduced infiltration. Conversely, lower coefficients in years like 2011 and 2020 suggest greater infiltration or retention of precipitation within the catchment. This variability highlights the influence of climatic and catchment conditions on runoff generation.

2.6. Calculate Ratio Between Evapotranspiration and Precipitation

# Group data by Year and calculate total evapotranspiration and total precipitation
evapotrans_prec <- main_data %>% 
  group_by(Year) %>% 
  summarise(
    Avg_Evapotrans = sum(Evapotrans),
    Avg_Precipitation = sum(Avg_Precipitation))

# Calculate evapotranspiration to precipitation ratio
evapotrans_prec$Ratio <- evapotrans_prec$Avg_Evapotrans / evapotrans_prec$Avg_Precipitation

evapotrans_prec

## # A tibble: 13 × 4
##     Year Avg_Evapotrans Avg_Precipitation Ratio
##    <dbl>          <dbl>             <dbl> <dbl>
##  1  2010           789.             1603. 0.492
##  2  2011           862.             1168. 0.739
##  3  2012           887.             1514. 0.586
##  4  2013           819.             1484. 0.552
##  5  2014           784.             1892. 0.414
##  6  2015           809.             1205. 0.672
##  7  2016           781.             1435. 0.544
##  8  2017           887.             1514. 0.586
##  9  2018           807.             1314. 0.614
## 10  2019           811.             1512. 0.537
## 11  2020           832.             1445. 0.576
## 12  2021           819.             1403. 0.584
## 13  2022           862.             1151. 0.749

# Plot evapotranpitation to precipitation ratio
plot(evapotrans_prec$Year, evapotrans_prec$Ratio, type = "b", col = "orange",
     xlab = "Year", ylab = "Ratio", main = "Evapotranspiration to Precipitation Ratio")

The evapotranspiration to precipitation ratio demonstrates annual variations, ranging from approximately 0.41 to 0.75. Years like 2011 and 2022 show higher ratios, indicating a larger fraction of precipitation is lost to evapotranspiration. Conversely, years like 2014 and 2010 have lower ratios, suggesting more precipitation is available for runoff and infiltration. These fluctuations highlight the influence of climatic variability on water balance dynamics within the catchment.

3. LOW FLOW ANALYSIS AND TREND ANALYSIS

This section focuses on low-flow and trend analysis of the hydrological data. We will conduct a low-flow analysis using the lfstat package, which includes baseflow separation and plotting a flow duration curve. Additionally, we will calculate key low-flow indexes such as BFI (Baseflow Index), MAM (Mean Annual Minimum flow), and quantile flows (Q95, Q90, Q70). The section will also explore drought periods based on the Q95 threshold and analyze seasonal trends, including the Mann-Kendall trend test for variables like discharge, precipitation, air temperature, and evapotranspiration.

3.1. Conduct Low-Flow Analysis

# Convert the Date column to POSIXct format for time series handling
main_data$Date <- as.POSIXct(strptime(main_data$Date, format = "%Y-%m-%d"))

# Define a zoo object for Discharge data with Date as the index
Qzoo <- zoo(main_data$Discharge, main_data$Date)

# Convert zoo object to lfobj using createlfobj
discharge_timeseries <- createlfobj(ts(Qzoo), startdate = "01/01/2010", hyearstart = 1)

# Define the units for the lfobj data
setlfunit("m^3/s")

# Conduct baseflow separation
discharge_timeseries$baseflow <- baseflow(discharge_timeseries$flow)

# Display the first few rows of the discharge_timeseries object
head(discharge_timeseries)

##   day month year   flow hyear baseflow
## 1   1     1 2010 59.058  2010       NA
## 2   2     1 2010 53.925  2010       NA
## 3   3     1 2010 47.681  2010       NA
## 4   4     1 2010 41.040  2010       NA
## 5   5     1 2010 38.694  2010       NA
## 6   6     1 2010 36.966  2010       NA

plot(discharge_timeseries$flow, type = "l", ylab = "Flow")
lines(discharge_timeseries$baseflow, col = "green")

The low-flow analysis plot displays the discharge over time (black line) along with the baseflow (green line) derived from the baseflow separation process. A noticeable decreasing trend in total flow is evident over the years, which may indicate long-term changes in climatic conditions, reduced precipitation, or increased water abstraction. The baseflow component remains relatively stable, indicating the consistent contribution of groundwater to the streamflow, even as overall flow declines.

bfplot(discharge_timeseries)

## NULL

3.2. Plot Flow Duration Curve

# Plot the Flow Duration Curve
fdc(discharge_timeseries)

##                  100%      99%      98%      97%      96%      95%      94%
## FDC-quantiles 859.744 265.1749 198.5197 165.5311 139.5318 124.5457 111.4355
##                    93%      92%      91%     90%     89%      88%      87%
## FDC-quantiles 101.1971 94.24816 87.01376 81.7343 75.5104 71.20372 67.21979
##                   86%      85%      84%      83%     82%      81%    80%
## FDC-quantiles 63.0349 59.76245 57.06852 54.55028 52.1816 50.30966 47.915
##                    79%      78%      77%     76%     75%      74%      73%
## FDC-quantiles 46.00891 44.48692 42.86936 41.7912 40.2515 39.03082 37.94931
##                    72%      71%     70%      69%      68%      67%      66%
## FDC-quantiles 36.84004 36.01759 35.1008 34.39782 33.47836 32.65884 31.92196
##                    65%      64%      63%      62%      61%    60%      59%
## FDC-quantiles 30.96445 30.20936 29.56432 28.97764 28.25367 27.558 26.99014
##                    58%      57%      56%     55%      54%      53%      52%
## FDC-quantiles 26.24132 25.62853 25.02784 24.5076 24.01942 23.52938 22.99832
##                    51%     50%      49%      48%      47%      46%      45%
## FDC-quantiles 22.40255 22.0315 21.56409 21.09428 20.65336 20.27288 19.79815
##                    44%      43%     42%      41%     40%      39%      38%
## FDC-quantiles 19.38068 18.96778 18.5832 18.25586 17.8788 17.58499 17.24476
##                    37%      36%      35%      34%      33%      32%     31%
## FDC-quantiles 16.95739 16.73312 16.48505 16.18994 15.93157 15.70604 15.4277
##                   30%      29%      28%      27%      26%      25%      24%
## FDC-quantiles 15.1942 14.95626 14.73848 14.46052 14.20366 13.96625 13.71984
##                    23%      22%      21%     20%      19%    18%    17%
## FDC-quantiles 13.48862 13.22834 12.97761 12.7024 12.50939 12.278 12.051
##                    16%      15%      14%      13%    12%      11%     10%
## FDC-quantiles 11.78856 11.53025 11.31492 11.05477 10.793 10.61953 10.3591
##                   9%      8%      7%      6%      5%      4%    3%      2%
## FDC-quantiles 10.027 9.69504 9.34229 8.95922 8.63805 8.29416 7.883 7.42382
##                    1%    0%
## FDC-quantiles 6.87029 5.515

abline(h=Q95(discharge_timeseries), col = "green", lty = 2)
text(90, 10, "Q95")

The Flow Duration Curve (FDC) shows the distribution of flow magnitudes over time in the catchment. The Q95 threshold (green dashed line) represents the flow that is exceeded 95% of the time, indicating low-flow conditions. The steep slope of the FDC at higher flows indicates a high variability of peak flows, while the flatter curve at lower flows reflects a stable baseflow regime. This result demonstrates the dominance of low to moderate flows in the catchment, with occasional high-flow events contributing to the variability observed in the upper range of the curve.

# FDC for the year 2014
y2014 <- fdc(discharge_timeseries, year = 2014)

# FDC for the year 2022
y2022 <- fdc(discharge_timeseries, year = 2022)

# Plot Flow Duration Curves (FDC) for 2014 and 2022
plot(as.numeric(y2014), type = "l", col = "blue", ylab = "Flow", main = "Flow Frequency Comparison: 2014 and 2022")
lines(as.numeric(y2022), col = "red", type = "l", lty = 2)
legend("topright", legend = c("2014", "2022"), col = c("blue", "red"), lty = c(1, 2))

The Flow Duration Curve (FDC) comparison between 2014 (blue line) and 2022 (red dashed line) shows a clear decrease in flow magnitudes over time. The FDC for 2014 indicates higher flows across all exceedance frequencies, suggesting wetter conditions and greater runoff compared to 2022. The lower flows in 2022 reflect a drier catchment, potentially due to reduced precipitation, increased water abstraction, or changes in land use. This comparison shows the progressive drying of the catchment and the variability in hydrological behavior over time.

3.3. Calculate Hydrological Indexes

# Calculate Baseflow Index
bfi <- BFI(discharge_timeseries)
bfi

## [1] 0.4925945

# Calculate Mean Annual Minimum Flow
mam <- MAM(discharge_timeseries)
mam

## [1] 8.581835

# Calculate Mean Flow
meanflow <- mean(discharge_timeseries$flow, na.rm = TRUE)
meanflow

## [1] 38.48014

# Calculate Flow Quantiles
# Q95
q95 <- Q95(discharge_timeseries)
q95

##     Q95 
## 8.63805

# Q90
q90 <- Q90(discharge_timeseries)
q90

##     Q90 
## 10.3591

# Q70
q70 <- Q70(discharge_timeseries)
q70

##     Q70 
## 15.1942

# Calculate Seasonality Index at Q95
seasonal_index <- seasindex(discharge_timeseries, Q = 95)
seasonal_index

## $theta
## [1] 3.857201
## 
## $D
## [1] 224.0708
## 
## $r
## [1] 0.7469678

# Calculate Seasonality Ratio
seasonal_ratio <- seasratio(discharge_timeseries, breakdays = "01/07", Q = 95)
seasonal_ratio

## Q95(1/1-01/07)/Q95(01/07-1/1) 
##                        1.3899

# Summarize 
summary(discharge_timeseries)

## Startdate:  2010-01-01     (calendar year)
## Enddate:    2022-12-31     (calendar year)
## 
## Time series covers 13 years.
## The hydrological year is set to start on January 1st.
## 
## Meanflow     MAM7      Q95      BFI 
##  38.4800   8.5820   8.6380   0.4926

3.4. Identify Drought Periods Based on Q95 Threshold

# Identify Drought Periods
discharge_timeseries$drought <- discharge_timeseries$flow < q95

# Display few rows of identified drought periods
drought_periods = discharge_timeseries[discharge_timeseries$drought == TRUE, ]
head(drought_periods)

##     day month year  flow hyear baseflow drought
## 191  10     7 2010 8.472  2010   8.4490    TRUE
## 192  11     7 2010 8.172  2010   8.1720    TRUE
## 193  12     7 2010 7.909  2010   7.9090    TRUE
## 194  13     7 2010 8.087  2010   7.8908    TRUE
## 195  14     7 2010 8.161  2010   7.8726    TRUE
## 196  15     7 2010 7.961  2010   7.8544    TRUE

# Create a Date column
drought_periods$Date <- as.Date(with(drought_periods, paste(year, month, day, sep = "-")), "%Y-%m-%d")

plot(drought_periods$Date, drought_periods$flow,
     col = "#1a24a4", xlab = "Date", ylab = "Flow", main = "Flow with Drought Periods (Below Q95 Threshold)")

The drought periods are identified based on the Q95 threshold (8.638 m³/s), representing flows below which drought conditions occur. The plot shows variability in drought occurrences across the years. Drought events are minimal in 2010 and 2011 but increase significantly in 2012 and 2013, with prolonged low-flow conditions. No droughts are observed between 2014 and 2016 or 2019 and 2021, indicating periods of higher flows. In contrast, 2017 and 2022 show moderate to high numbers of drought days, reflecting a return to frequent low-flow conditions.

# Identify drought periods using the Q95 threshold
drought_time = find_droughts(discharge_timeseries, threshold = Q95(discharge_timeseries), na.rm = TRUE)

## Warning in as.xts.lfobj(x): No unit found in attributes, assuming 'm³/s'.
## Use flowunit(x) <- "l/s" to define the flow unit. See help(flowunit).

# Summarize drought events, filtering out minor ones
selected = summary(drought_time, drop_minor = 0)

# Plot the drought time series
plot(drought_time)

The plot displays drought periods identified using the Q95 threshold, with shaded segments indicating periods of low flow. The results show shorter drought durations around the end of 2010, followed by multiple droughts between 2012 and 2014, some of which are more prolonged. A significant drought period is observed toward the end of 2017, with even longer durations in the year 2022.

# Plot the seasonal distribution of selected drought events
plot(season(selected$time))

The plot illustrates the seasonal distribution of drought events. The results indicate that droughts predominantly occur during the summer season, followed by autumn. Spring experiences fewer droughts, while winter has the least occurrences due to lower evapotranspiration and stable baseflow conditions during colder months.

3.5. Conduct Mann-Kendall Test

# Initialize list to store the results
mann_kandell_test <- list()

# Mann-Kendall test for Discharge
mann_kandell_test$Discharge <- MannKendall(main_data$Discharge)

# Mann-Kendall test for Precipitation at each station
precip_stations <- c("P1_Luce", "P2_Solcava", "P3_Lasko", "P4_GornjiGrad", "P5_Radegunda")
for (station in precip_stations) {
  mann_kandell_test[[station]] <- MannKendall(main_data[[station]])
}

# Mann-Kendall test for Air Temperature
mann_kandell_test$Airtemp <- MannKendall(main_data$Airtemp)

# Mann-Kendall test for Evapotranspiration
mann_kandell_test$Evapotrans <- MannKendall(main_data$Evapotrans)

# Test result
mann_kandell_test

## $Discharge
## tau = -0.0314, 2-sided pvalue =0.0011597
## 
## $P1_Luce
## tau = -0.0207, 2-sided pvalue =0.051273
## 
## $P2_Solcava
## tau = -0.00822, 2-sided pvalue =0.44307
## 
## $P3_Lasko
## tau = 0.000218, 2-sided pvalue =0.98394
## 
## $P4_GornjiGrad
## tau = -0.0113, 2-sided pvalue =0.28991
## 
## $P5_Radegunda
## tau = -0.00613, 2-sided pvalue =0.56852
## 
## $Airtemp
## tau = 0.0207, 2-sided pvalue =0.032784
## 
## $Evapotrans
## tau = 0.00605, 2-sided pvalue =0.53637

The Mann-Kendall test results indicate no significant trends for most precipitation stations, with P1_Luce showing a marginally decreasing trend. For air temperature, a slight increasing trend is observed with a statistically significant p-value, suggesting a gradual warming over time. Evapotranspiration shows no significant trend, indicating stable evapotranspiration levels throughout the analyzed period.

mk.test(main_data$Discharge)

## 
##  Mann-Kendall trend test
## 
## data:  main_data$Discharge
## z = -3.2486, n = 4748, p-value = 0.00116
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
## -3.543340e+05  1.189671e+10 -3.144278e-02

The Mann-Kendall test for discharge indicates a statistically significant decreasing trend, with a p-value of 0.00116. This result reflects a decline in water flow over time, suggesting potential long-term changes in hydrological conditions or water availability in the catchment.

plot(main_data$Date, main_data$Discharge, type = "l", xlab = "Year", ylab = "Discharge[m3/s]", main = "Discharge plot")
abline(lm(main_data$Discharge ~ main_data$Date), lty = 2, col = "red")

4. FLOOD FREQUENCY ANALYSIS

In this section, we will perform a flood frequency analysis based on daily discharge data to assess extreme flood events. The analysis begins with defining an annual maxima sample derived from daily discharge records. Next, we will import the vQvk peak discharge data, covering the complete available dataset, and compare it with the annual maxima sample to validate consistency. Using the vQvk data, we will conduct a flood frequency analysis applying at least three different distribution functions, and the results will be visualized, highlighting measured discharge data under the Weibull distribution. Additionally, confidence intervals will be calculated by the genci.simple function from the lmomco package.

4.1. Define the Annual Maxima Sample

# Max discharge values 
DischargeMax <- main_data %>% 
  group_by(Year) %>% 
  summarise(
    maxDischarge = max(Discharge)
  )

DischargeMax

## # A tibble: 13 × 2
##     Year maxDischarge
##    <dbl>        <dbl>
##  1  2010         860.
##  2  2011         323.
##  3  2012         579.
##  4  2013         443.
##  5  2014         571.
##  6  2015         397.
##  7  2016         402.
##  8  2017         582.
##  9  2018         282.
## 10  2019         340.
## 11  2020         322.
## 12  2021         243.
## 13  2022         330.

# Plot yearly maximum discharge
plot(DischargeMax$Year, DischargeMax$maxDischarge, type = "l", xlab = "Year", ylab = "Max Discharge[m3/s]", 
     main = "Yearly Max Discharge")
abline(lm(DischargeMax$maxDischarge ~ DischargeMax$Year), lty = 2, col = "red")

The plot of yearly maximum discharge shows a clear declining trend in peak discharge values over the analyzed period. This is supported by the red dashed line, which represents the linear regression fit, further confirming the decreasing trend. The reduction in maximum discharge values could indicate a potential change in hydrological conditions, such as reduced precipitation intensity or changes in land use and catchment characteristics, affecting runoff generation during extreme events.

4.2. Import the vQvk Peak Discharge Data and Compare Values

peak_data <- read_excel("peak_discharge_data.xlsx")

# Calculate annual peak for vQvk discharge
annual_maxima_vQvk <- peak_data %>% 
  group_by(Year) %>% 
  summarise(
    vQvk = max(Qvk, na.rm = TRUE)
  )

# Display few rows
head(annual_maxima_vQvk)

## # A tibble: 6 × 2
##   Year   vQvk
##   <chr> <dbl>
## 1 1955    421
## 2 1956    569
## 3 1957    290
## 4 1958    524
## 5 1959    484
## 6 1960    502

# Filter vQvk data from 2010 onward
annual_vQvk_selected <- annual_maxima_vQvk %>% 
  filter(Year >= 2010)

# Combine the two peak discharge
combined_data <- data.frame(
  Year = annual_vQvk_selected$Year,
  main_DischargeMax = DischargeMax$maxDischarge,
  vQvk_Discharge = annual_vQvk_selected$vQvk
)

# Determine y-axis limits
y_range <- range(c(combined_data$main_DischargeMax, combined_data$vQvk_Discharge), na.rm = TRUE)

# Plotting vQvk peak discharge
plot(combined_data$Year, combined_data$vQvk_Discharge, type = "l", col = "blue", lwd = 1,
     xlab = "Year", ylab = "Peak Discharge", main = "Comparison of Peak Discharge Values", ylim = y_range)

# Add main_DischargeMax to the plot
lines(combined_data$Year, combined_data$main_DischargeMax, col = "red", lwd = 1)

# Adding a legend
legend("topright", legend = c("Main DischargeMax", "vQvk Discharge"),
       col = c("red", "blue"), lwd = 1)

4.3. Conduct Flood Frequency Analysis Using Multiple Distribution Functions

# Summary statistics for vQvk
summary(annual_maxima_vQvk)

##      Year                vQvk       
##  Length:66          Min.   : 262.3  
##  Class :character   1st Qu.: 461.7  
##  Mode  :character   Median : 565.0  
##                     Mean   : 644.3  
##                     3rd Qu.: 766.8  
##                     Max.   :1490.0

# Plot annual maxima series for vQvk
plot(annual_maxima_vQvk$Year, annual_maxima_vQvk$vQvk, type = "b",
     xlab = "Year", ylab = "vQvk [m3/s]", main = "Annual Peak Discharge for Savinja River")

# Calculate L-moments for the vQvk data
lmoments <- lmoms(annual_maxima_vQvk$vQvk)

# Display calculated L-moments
lmoments$lambdas  

## [1] 644.270955 146.790664  41.173791  27.376282   5.118641

# Estimate parameters for PE3, GEV, GNO, Normal, and Weibull distributions
pe3par <- lmom2par(lmoments, type = "pe3")
gevpar <- lmom2par(lmoments, type = "gev")
gnopar <- lmom2par(lmoments, type = "gno")
norpar <- lmom2par(lmoments, type = "nor")  
wei_par <- lmom2par(lmoments, type = "wei") 

# Define return periods and calculate exceedance probabilities
return_periods <- c(2, 5, 10, 20, 30, 50, 100, 300,500)

# Calculate Fx values for each return period
exceedance_probs <- 1 - 1 / return_periods  

# Calculate discharge values for each return period using PE3, GEV, GNO, Normal, and Weibull distributions
discharge_pe3 <- quape3(exceedance_probs, pe3par)
discharge_gev <- quagev(exceedance_probs, gevpar)
discharge_gno <- quagno(exceedance_probs, gnopar)
discharge_nor <- quanor(exceedance_probs, norpar)
discharge_wei <- quawei(exceedance_probs, wei_par)

# Calculate Weibull plotting positions for observed data and return periods for measured data
weibull_positions <- pp(annual_maxima_vQvk$vQvk, a = 0)
observed_return_periods <- 1 / (1 - weibull_positions)

# Plot the PE3, GEV, GNO, Normal, and Weibull distribution results
plot(return_periods, discharge_pe3, type = "l", log = "x", col = "blue",
     xlab = "Return Period [years]", ylab = "Discharge [m3/s]",
     main = "Flood Frequency Analysis for Savinja River", ylim = c(400, max(discharge_pe3, discharge_gev, discharge_gno, discharge_wei, na.rm = TRUE)))

# Add GEV, GNO, Normal, and Weibull distribution lines
lines(return_periods, discharge_gev, col = "red")
lines(return_periods, discharge_gno, col = "green")
lines(return_periods, discharge_nor, col = "orange")
lines(return_periods, discharge_wei, col = "#ff31a1") 

# Add measured discharge data points
points(observed_return_periods, sort(annual_maxima_vQvk$vQvk), col = "black", pch = 16)

# Add legend
legend("topleft", legend = c("PE3 Distribution", "GEV Distribution", "GNO Distribution", 
                                 "Normal Distribution", "Weibull Distribution", "Observed Data"),
       col = c("blue", "red", "green", "orange", "#ff31a1", "black"), lty = 1, pch = c(NA, NA, NA, NA, NA))

The flood frequency analysis illustrates discharge values for various return periods using multiple probability distributions, including PE3, GEV, GNO, Normal, and Weibull. The Weibull and PE3 distributions show a closer fit to the observed data, particularly for low and moderate return periods, making them more reliable for this dataset. On the other hand, GEV and GNO distributions tend to overestimate discharge values, especially for higher return periods, while the Normal distribution consistently underestimates discharge across most return periods. This comparison shows the importance of selecting an appropriate distribution for accurate flood frequency analysis in hydrological studies.

4.4. Calculate Confidence Intervals

# Set seed for reproducibility
set.seed(123)

# Calculate confidence intervals for the PE3 distribution
ci_pe3 <- genci.simple(para = pe3par, n = length(annual_maxima_vQvk$vQvk), 
                       f = exceedance_probs, level = 0.90, nsim = 1000, edist = "gno")

## 9-8-7-6-5-4-3-2-1-0-

# Combine the confidence intervals and design discharges
confidence_intervals <- data.frame(
  Return_Period = return_periods,
  Lower_CI = ci_pe3$lwr,
  Upper_CI = ci_pe3$upr,
  PE3_Discharge = discharge_pe3
)

# Plot PE3 distribution with confidence intervals
plot(confidence_intervals$Return_Period, confidence_intervals$PE3_Discharge, type = "l", log = "x", 
     col = "green", xlab = "Return Period [years]", ylab = "Discharge [m3/s]",
     main = "Flood Frequency Analysis", 
     ylim = c(400, max(ci_pe3$upr, na.rm = TRUE)))

# Add confidence interval bounds for PE3 distribution
lines(confidence_intervals$Return_Period, confidence_intervals$Lower_CI, col = "purple", lty = 2)
lines(confidence_intervals$Return_Period, confidence_intervals$Upper_CI, col = "purple", lty = 2)

# Add observed discharge data points for comparison
points(observed_return_periods, sort(annual_maxima_vQvk$vQvk), col = "black", pch = 16)

# Add legend
legend("bottomright", legend = c("PE3 Distribution", "Lower CI", "Upper CI", "Observed Data"),
       col = c("green", "purple", "purple", "black"), lty = c(1, 2, 2, NA), pch = c(NA, NA, NA, 16))

Based on the previous analysis, the Weibull and PE3 distributions were observed to fit the dataset well, with PE3 being slightly more consistent with the observed data. As a result, the PE3 distribution was selected for calculating the confidence intervals. The plot shows the PE3 discharge curve (green line) along with the 90% confidence intervals (purple dashed lines). The observed data points (black dots) are mostly contained within the confidence intervals, demonstrating the reliability of the PE3 model. Additionally, the confidence intervals widen for larger return periods, reflecting increased uncertainty in predicting extreme events such as rare floods.

6. PRECIPITATION DATA ANALYSIS

In this section, we will analyze precipitation data to explore trends and variability across multiple thresholds and indices. We begin by calculating the total number of days with measurable precipitation, as well as days with precipitation exceeding 10 mm and 50 mm, followed by determining the probability of occurrence for each threshold. Next, we will compute Standardized Precipitation Index (SPI) values over 1, 3, 6, and 12-month periods for each of the five precipitation stations to assess drought conditions. The identified drought periods will then be compared to those derived from discharge data using the lfstat package, providing insight into the relationship between precipitation and drought events. Additionally, we will calculate the Effective Drought (ED) index and examine its relationship with the SPI results. Finally, we will incorporate the North Atlantic Oscillation (NAO) index by downloading NAO data and calculating its correlation with precipitation records to investigate potential climate influences on regional precipitation patterns.

6.1 Calculate Days with Precipitation

# Calculate days with precipitation > 10 mm
prep_above_10 <- main_data %>%
  summarise(
    P1_Luce = sum(P1_Luce > 10),
    P2_Solcava = sum(P2_Solcava > 10),
    P3_Lasko = sum(P3_Lasko > 10),
    P4_GornjiGrad = sum(P4_GornjiGrad > 10),
    P5_Radegunda = sum(P5_Radegunda > 10),
  )

# Calculate frequency for days with precipitation > 10 mm
frequency_above_10 <- prep_above_10 * 100 / nrow(main_data)

# Combine both the counts and frequency rows
results_above_10 <- bind_rows(prep_above_10, frequency_above_10)
rownames(results_above_10) <- c("Days", "Frequency")

# Display the results
results_above_10

##             P1_Luce P2_Solcava  P3_Lasko P4_GornjiGrad P5_Radegunda
## Days      624.00000  627.00000 538.00000     640.00000    589.00000
## Frequency  13.14238   13.20556  11.33109      13.47936     12.40522

# Calculate days with precipitation > 50 mm
prep_above_50 <- main_data %>%
  summarise(
    P1_Luce = sum(P1_Luce > 50),
    P2_Solcava = sum(P2_Solcava > 50),
    P3_Lasko = sum(P3_Lasko > 50),
    P4_GornjiGrad = sum(P4_GornjiGrad > 50),
    P5_Radegunda = sum(P5_Radegunda > 50)
  )

# Calculate frequency for days with precipitation > 50 mm
frequency_above_50 <- prep_above_50 * 100 / nrow(main_data)

# Combine both the counts and frequency rows
results_above_50 <- bind_rows(prep_above_50, frequency_above_50)
rownames(results_above_50) <- c("Days", "Frequency")

# Display the results
results_above_50

##             P1_Luce P2_Solcava  P3_Lasko P4_GornjiGrad P5_Radegunda
## Days      59.000000  62.000000 18.000000     50.000000   31.0000000
## Frequency  1.242628   1.305813  0.379107      1.053075    0.6529065

6.2 Calculate SPI Indexes

# Create a zoo object for all 5 stations
precip_zoo <- zoo(main_data[, c("P1_Luce", "P2_Solcava", "P3_Lasko", "P4_GornjiGrad", "P5_Radegunda")], order.by = main_data$Date)

head(precip_zoo)

##            P1_Luce P2_Solcava P3_Lasko P4_GornjiGrad P5_Radegunda
## 2010-01-01     7.0        4.8      4.3           4.6          2.2
## 2010-01-02     3.0        0.6      2.1           3.2          2.0
## 2010-01-03     0.1        2.7      1.5           0.0          0.1
## 2010-01-04     0.0        0.0      0.0           0.0          0.0
## 2010-01-05     3.3        2.0      7.0           3.1          4.0
## 2010-01-06     6.2        4.6      7.0           9.2          8.2

# Convert time index to year-month
yearmon_index <- as.yearmon(time(precip_zoo))

# Aggregate precipitation data by month
monthly_precip <- aggregate(precip_zoo, yearmon_index, sum)

# Plot the aggregated monthly precipitation data for P1_Luce
plot(monthly_precip$P1_Luce, xlab="Year", ylab="Monthly precipitation [mm]", main="Monthly Precipitation P1_Luce")

# Get summary of the monthly data
summary(monthly_precip)  

##      Index         P1_Luce        P2_Solcava       P3_Lasko     
##  Min.   :2010   Min.   :  0.1   Min.   :  0.3   Min.   :  0.00  
##  1st Qu.:2013   1st Qu.: 66.5   1st Qu.: 61.0   1st Qu.: 55.90  
##  Median :2016   Median :117.3   Median :112.4   Median : 87.50  
##  Mean   :2016   Mean   :128.1   Mean   :128.0   Mean   : 98.45  
##  3rd Qu.:2020   3rd Qu.:174.7   3rd Qu.:185.1   3rd Qu.:135.80  
##  Max.   :2023   Max.   :382.7   Max.   :398.0   Max.   :313.20  
##  P4_GornjiGrad    P5_Radegunda  
##  Min.   :  0.0   Min.   :  0.0  
##  1st Qu.: 67.2   1st Qu.: 63.6  
##  Median :113.4   Median :103.0  
##  Mean   :124.6   Mean   :114.5  
##  3rd Qu.:166.6   3rd Qu.:166.4  
##  Max.   :360.6   Max.   :326.0

# Calculate SPI at 3-month scale for all 5 stations
# P1_Luce SPI (3-month scale)
spi_3_P1_Luce <- spi(data = as.numeric(monthly_precip$P1_Luce), scale = 3, distribution = "Gamma")

## [1] "Calculating the Standardized Precipitation Evapotranspiration Index (SPEI) at a time scale of 3. Using kernel type 'rectangular', with 0 shift. Fitting the data to a Gamma distribution. Using the ub-pwm parameter fitting method. Checking for missing values (`NA`): all the data must be complete. Using the whole time series as reference period. Input type is vector. No time information provided, assuming a monthly time series."

# P2_Solcava SPI (3-month scale)
spi_3_P2_Solcava <- spi(data = as.numeric(monthly_precip$P2_Solcava), scale = 3, distribution = "Gamma")

## [1] "Calculating the Standardized Precipitation Evapotranspiration Index (SPEI) at a time scale of 3. Using kernel type 'rectangular', with 0 shift. Fitting the data to a Gamma distribution. Using the ub-pwm parameter fitting method. Checking for missing values (`NA`): all the data must be complete. Using the whole time series as reference period. Input type is vector. No time information provided, assuming a monthly time series."

# P3_Lasko SPI (3-month scale)
spi_3_P3_Lasko <- spi(data = as.numeric(monthly_precip$P3_Lasko), scale = 3, distribution = "Gamma")

## [1] "Calculating the Standardized Precipitation Evapotranspiration Index (SPEI) at a time scale of 3. Using kernel type 'rectangular', with 0 shift. Fitting the data to a Gamma distribution. Using the ub-pwm parameter fitting method. Checking for missing values (`NA`): all the data must be complete. Using the whole time series as reference period. Input type is vector. No time information provided, assuming a monthly time series."

# P4_GornjiGrad SPI (3-month scale)
spi_3_P4_GornjiGrad <- spi(data = as.numeric(monthly_precip$P4_GornjiGrad), scale = 3, distribution = "Gamma")

## [1] "Calculating the Standardized Precipitation Evapotranspiration Index (SPEI) at a time scale of 3. Using kernel type 'rectangular', with 0 shift. Fitting the data to a Gamma distribution. Using the ub-pwm parameter fitting method. Checking for missing values (`NA`): all the data must be complete. Using the whole time series as reference period. Input type is vector. No time information provided, assuming a monthly time series."

# P5_Radegunda SPI (3-month scale)
spi_3_P5_Radegunda <- spi(data = as.numeric(monthly_precip$P5_Radegunda), scale = 3, distribution = "Gamma")

## [1] "Calculating the Standardized Precipitation Evapotranspiration Index (SPEI) at a time scale of 3. Using kernel type 'rectangular', with 0 shift. Fitting the data to a Gamma distribution. Using the ub-pwm parameter fitting method. Checking for missing values (`NA`): all the data must be complete. Using the whole time series as reference period. Input type is vector. No time information provided, assuming a monthly time series."

# Plot SPI-3 for P1_Luce
plot(spi_3_P1_Luce, main = "SPI-3 for P1_Luce")

## Warning: Removed 155 rows containing missing values (`geom_point()`).

# Plot SPI-3 for P2_Solcava
plot(spi_3_P2_Solcava, main = "SPI-3 for P2_Solcava")

## Warning: Removed 155 rows containing missing values (`geom_point()`).

# Plot SPI-3 for P3_Lasko
plot(spi_3_P3_Lasko, main = "SPI-3 for P3_Lasko")

## Warning: Removed 155 rows containing missing values (`geom_point()`).

# Plot SPI-3 for P4_GornjiGrad
plot(spi_3_P4_GornjiGrad, main = "SPI-3 for P4_GornjiGrad")

## Warning: Removed 155 rows containing missing values (`geom_point()`).

# Plot SPI-3 for P5_Radegunda
plot(spi_3_P5_Radegunda, main = "SPI-3 for P5_Radegunda")

## Warning: Removed 155 rows containing missing values (`geom_point()`).

# View SPI values
spi_3_P1_Luce

##   [1]            NA            NA -0.7204634761 -0.1035448200 -0.5539357550
##   [6] -0.6676236843 -1.2969851317 -0.0655589479  0.3259998529  1.7653594567
##  [11]  1.1743821379  1.0562537491  0.7864004975  0.5315218041  0.0836068510
##  [16] -0.3778538826 -0.1503402766  0.0434250727  0.0874669662  0.5337545988
##  [21]  0.4556579062 -0.4739078884  0.1009485414 -1.0554242865 -0.5487552819
##  [26] -1.9509528509 -0.7768173785 -1.8320052570 -1.2106029338  0.0593617193
##  [31]  0.7865403272  0.0111125757  0.4313171070  1.0155799377  1.8984148098
##  [36]  2.0180754590  1.6783405425  0.3900624293  0.2186045794  0.9304389929
##  [41]  1.1681300011  1.5943340047 -0.2344968859 -0.7152876314 -2.0219231681
##  [46] -1.6088188273 -0.8308576836  0.2313359447  0.3006416530  1.4541264415
##  [51]  1.8924172585  1.7712864233  1.3287712361 -0.9517816819 -0.1692218462
##  [56]  0.4717715631  1.7397842350  1.5117932583  0.7086622877  0.8132757906
##  [61]  0.5325371731  0.5917726450 -0.6862657627 -0.1845212838 -1.0187620935
##  [66] -0.5429304619 -0.2155548230  1.1056130145  0.4380425893  0.0857337977
##  [71] -0.0005738679 -0.5370730805 -2.0228457110 -1.2656144437  0.7686494723
##  [76]  1.2674089475  1.0134699039 -0.6940074221 -0.6938430762 -0.0260690707
##  [81]  0.6011259267 -0.7656586966 -0.8326843161 -0.1094490102 -0.0634693010
##  [86] -0.2205900906 -0.8041570121 -0.3723307295  1.1224227183  0.0385131016
##  [91]  1.0044809018 -1.7667037986 -0.6991236156  0.3616881567  0.1583698010
##  [96]  0.7313631578  0.8082305226  0.9764092742  0.9127374322  0.5074158067
## [101]  0.7981022991  1.0834361643  0.7710131441 -0.0281412397  0.0243632898
## [106] -0.6675149485  0.2305222973 -0.5486176001 -1.1294044301 -1.4507829219
## [111] -0.6551449107 -0.2509881927  0.0996966452  0.3182928461  0.1555393112
## [116]  0.1787635734 -0.2419443280 -0.2190713667 -1.0221250319  0.1876205391
## [121]  0.7997899101  0.6987957449 -0.2889669808 -0.3224278584 -0.6310426306
## [126] -0.4581669440 -0.0433886664  0.7751208418  1.0284666199  0.0982773229
## [131]  0.5755299204 -0.7070806847  0.3580075641  0.4135041528  1.2872524581
## [136]  0.4741089368 -0.6088379641  1.6954463058  1.8524733713  1.6152503773
## [141] -0.4130392476 -0.6971142888 -1.5138998942 -1.2742193808 -0.6216554028
## [146] -0.1473223838 -0.7668758330 -1.2083643412 -1.0946375692 -1.3541367833
## [151] -1.9365516846 -2.0665831884 -1.6110726896 -0.2567218298 -0.4366623883
## [156] -0.6084039781 -0.7632183398

spi_3_P2_Solcava

##   [1]           NA           NA -0.709959804 -0.174533764 -0.497014226
##   [6] -0.430933635  0.011703359 -0.002330165  0.482201164  0.897560815
##  [11]  1.009847864  0.855183314  0.678782311  0.335845732 -0.121475490
##  [16] -0.440607366 -0.323648029  0.215482312  0.756181542  1.017861769
##  [21]  0.507424063 -0.170113402  0.073484103 -0.986796098 -1.101574997
##  [26] -1.685578121 -0.840649208 -1.767250300 -1.317421508 -0.280943819
##  [31]  0.126513193 -0.238080607  0.149988713  0.986911951  1.645075609
##  [36]  1.783900628  1.110181653  0.172371864  0.028718380  0.805148911
##  [41]  0.928405189  1.732601397  0.030641832 -0.198973539 -1.700147450
##  [46] -1.034814771 -0.666018796  0.580198655  0.757928699  1.443396110
##  [51]  1.881046606  1.968251405  1.358931228 -0.333713907  1.044196045
##  [56]  1.090161763  1.702052692  1.263730411  0.692481686  0.985183101
##  [61]  0.388022824  0.574388943 -0.416976499  0.132888314 -0.609506906
##  [66]  0.035820221  0.344004865  1.160867857  0.251307581  0.145834550
##  [71]  0.264800164 -0.276429035 -1.785874198 -1.499691736  0.385879799
##  [76]  1.057238500  1.057162577 -0.528289628 -0.771281198 -0.760759640
##  [81] -0.095289536 -1.464799964 -1.152047817 -0.817679789 -0.688842188
##  [86] -0.467649344 -1.023298834 -0.584115626  0.895222854 -0.355491530
##  [91]  0.216427568 -2.139273378 -0.727013158  0.413397043  0.173568293
##  [96]  0.575557501  0.833935448  0.918890410  0.980835666  0.327726682
## [101]  0.496487980  0.467792171 -0.200755622 -0.684921188 -0.292476496
## [106] -1.181152146 -0.072779909 -0.955516915 -1.005298359 -1.203718306
## [111] -0.194040103  0.234508478  1.031230752  1.618783720  1.495363490
## [116]  0.312359304 -0.842060763 -0.577345665 -1.165186252  0.709291848
## [121]  1.214961313  0.960651077 -0.273540815 -0.340352115 -0.684128580
## [126] -1.096569331 -0.957714796  0.863493787  1.727339804  1.220022235
## [131]  1.170118347 -0.278758781  0.872198149  0.633472779  1.471187346
## [136]  0.330420936 -0.835875828  0.753879399  0.686983065  0.778310548
## [141] -0.555754956 -0.818326861 -1.264092922 -1.237272905 -0.510860093
## [146] -0.052460659 -0.604088821 -1.230098127 -1.273425239 -1.668011192
## [151] -2.823708556 -1.119114300 -0.462859516  0.495853942 -0.476376599
## [156] -0.743817312 -0.644799247

spi_3_P3_Lasko

##   [1]           NA           NA -0.318730311  0.488293739  0.065516561
##   [6]  0.140802544 -0.266518721 -0.148910456 -0.106848116  1.409528093
##  [11]  1.510970123  1.456046294  0.766348515  0.638035797 -0.852248834
##  [16] -0.829956754 -0.991936943 -0.828203542 -0.293253544  0.529311600
##  [21] -0.031549226 -1.359556061 -1.257381728 -1.281806457 -0.595221145
##  [26] -1.436736644 -0.521976870 -1.459844374 -0.725496543 -0.260814555
##  [31]  0.330452767 -0.598891737 -0.585915291  0.031033203  1.496063092
##  [36]  1.255116717  1.358655303  0.268482679  1.230483594  1.539408872
##  [41]  1.527435927  1.431104275 -0.575285249 -1.580625020 -2.284673982
##  [46] -1.047000350 -0.570143401  0.842604951  0.500459791  1.491209131
##  [51]  1.285053884  0.992755907  1.170389096 -0.120553861  1.158970685
##  [56]  0.573431040  1.373631969  1.044575667  0.558551282 -0.109421466
##  [61] -0.557813884  0.112793078  0.215496499  0.016989098 -0.703102756
##  [66] -0.283211706  1.088779137  0.955346370  0.355488130 -1.111079187
##  [71] -0.009017051 -0.249145716 -0.708948391 -1.433410298  0.831680932
##  [76]  1.209465025  1.306399922  0.650537236  0.791066622  0.086881572
##  [81]  0.207962502 -0.968355731 -0.613117465 -0.352262698 -0.362808855
##  [86] -0.346840636 -1.010045864 -0.346657549  0.078569339 -1.251390684
##  [91] -1.265874200 -2.300301324 -1.362828095  0.386930099  0.680867061
##  [96]  1.221848345  1.168069474  1.393903874  1.343450141  0.952756451
## [101]  0.967233236  0.507392742  0.409086977  0.210255956  0.305460454
## [106] -0.833081099 -1.337074546 -1.693878254 -2.922012442 -1.415180762
## [111] -1.280614179 -0.199270202 -0.129750387  1.176577925  1.127680993
## [116]  0.890921681  0.717808484  0.480683148 -0.002944508 -0.198057651
## [121]  0.322843367  0.219324308 -0.808543915 -1.022561747 -1.261416217
## [126] -1.715441140 -2.355398715  0.795597006  1.561060027  1.146513734
## [131]  0.516280430 -0.393488565  0.792705714  0.305967880  0.573673561
## [136] -0.204531595 -0.490676795  1.174080579  0.332098985  1.253281814
## [141]  0.320224062  0.306539246 -1.056993646 -0.591662718  0.138702313
## [146]  0.321065691 -0.386549797 -0.686029187 -0.457444061 -0.446384696
## [151] -0.447952062 -0.659201058 -0.432661621  0.680414206  0.278063602
## [156]  0.286271663  0.001287518

spi_3_P4_GornjiGrad

##   [1]           NA           NA -0.717472654  0.076834022 -0.509205904
##   [6] -0.641802564 -1.137559357 -0.277074732  0.052024124  1.730623873
##  [11]  1.401837588  1.465811714  0.918222709  0.629490871 -0.092858565
##  [16] -0.487177260 -0.772731759 -0.330128553 -0.910541763  0.312566750
##  [21] -0.085461014 -1.160102389 -0.434884080 -1.443658177 -0.660788358
##  [26] -1.731770298 -0.665023775 -1.612549391 -1.052019598 -0.537461031
##  [31]  0.100328925 -0.405067227  0.319167925  0.337925852  1.701775205
##  [36]  1.819079628  1.886701377  0.430722486  0.307797916  0.945729254
##  [41]  1.086669978  1.590199261 -0.419707846 -0.809905387 -1.668256578
##  [46] -1.188379664 -0.709739296  0.399384613  0.086265904  1.419926217
##  [51]  1.756156175  1.661769776  1.077692434 -1.173949001  0.350411446
##  [56]  0.775397181  2.031604855  1.492024916  0.441760628  0.420299822
##  [61]  0.170550693  0.561980492 -0.540429925 -0.384328244 -1.089042964
##  [66] -0.273592115 -0.120395372  0.846604922 -0.359493840 -0.516441045
##  [71]  0.017549446 -0.364089955 -1.546727260 -1.457216820  0.703219369
##  [76]  1.213283386  1.239926852  0.203678309  0.448897549  0.694729658
##  [81]  0.307810096 -1.101537535 -1.031462389 -0.293718858 -0.389226690
##  [86] -0.507335114 -0.767202848 -0.256958543  1.037415325 -0.039477339
##  [91]  1.597655885 -1.868749127 -0.459847253  0.456893742  0.189610133
##  [96]  0.646753503  0.785376096  1.010224670  1.154860798  0.746574675
## [101]  1.177492736  1.453076062  0.302829083  0.286010087  0.199736572
## [106] -0.232931152 -0.111420133 -0.614940656 -1.170389230 -1.260628318
## [111] -0.738937044 -0.252535653  0.084382412  0.373797320  1.076504481
## [116]  0.854066143  0.184266943 -0.046778143 -1.033464104  0.089868710
## [121]  0.679846076  0.655953201 -0.257905162 -0.247386060 -0.536325046
## [126] -0.541770764 -0.804365514  0.558123086  1.386605466  0.716869443
## [131]  1.042993738 -0.481366704  0.541156220  0.526610814  1.285417945
## [136]  0.427331563 -0.376722118  1.366523963  1.156068197  1.170155741
## [141] -0.154997473 -0.549505026 -1.247701950 -1.327269842 -0.340562236
## [146] -0.207353722 -1.014767545 -1.507194995 -1.150485306 -1.286264408
## [151] -1.540854110 -2.148708436 -1.699026634  0.233382711 -0.004262021
## [156] -0.131077374 -0.758694226

spi_3_P5_Radegunda

##   [1]          NA          NA -0.46053859  0.19529862 -0.56052002 -1.38783093
##   [7] -2.18017336 -1.39502840 -0.23133538  1.43361604  1.41272441  1.25066023
##  [13]  0.50004722  0.32821903 -0.34139795 -0.57117429 -0.66483260 -0.35792566
##  [19]  0.01956341  0.37962939  0.16036838 -1.18372706 -0.51010024 -1.35839118
##  [25] -0.71174862 -1.66791004 -0.77175665 -1.47716360 -1.12744349 -0.53879633
##  [31] -0.34788656 -0.40451088  0.36065990  0.62432914  1.77088452  1.62800829
##  [37]  1.79809306  0.47194107  0.65040389  1.38141828  1.64426005  1.65956725
##  [43] -0.80097369 -1.55638375 -2.36227460 -1.12382972 -0.37509596  0.94221759
##  [49]  0.76132537  1.60930134  1.77986025  1.42473820  1.02569029 -0.64434628
##  [55]  0.46654481  0.57965064  1.77813259  1.47034757  0.75836076  0.59761154
##  [61]  0.39920572  0.57927831 -0.52586587 -0.18738369 -0.61794857  0.12206524
##  [67]  0.51922133  1.02863770 -0.03497377 -0.79368756 -0.30577124 -0.44947984
##  [73] -1.36593051 -1.38925051  0.71747277  1.10143005  1.11323843 -0.29095495
##  [79]  0.57756668  0.69182793  0.54327157 -1.37012381 -1.15834451 -0.65424824
##  [85] -0.41275782 -0.69167089 -0.87962541 -0.58650300  0.62730506 -0.33372181
##  [91]  1.81379217 -0.58206517  0.02491121  0.33043076  0.26588199  0.83927758
##  [97]  0.46424837  0.88906381  1.04890927  0.86248803  1.15152911  1.27242044
## [103]  0.52072197 -0.37928760  0.11726215 -0.01987118  0.37430108 -0.83615355
## [109] -1.67353911 -1.22577520 -0.98860032 -0.23849151 -0.09396873  0.68324493
## [115]  0.34958355  0.83340617  0.12858089  0.26863449 -0.93984337  0.11954809
## [121]  0.81382205  0.70853189 -0.08534530 -0.35068384 -0.73831988 -0.83540773
## [127] -0.96464894  1.24462259  1.41200438  0.76549531  0.44217202 -0.45960842
## [133]  0.59028135  0.40284506  1.13664198  0.19383232 -0.65846001  1.44140306
## [139]  0.85236465  0.87667259 -0.92922779 -0.52770821 -1.33014400 -0.93796357
## [145] -0.19290923  0.09338470 -0.86253475 -1.33605039 -0.84822271 -0.62448379
## [151] -0.78182356 -1.22239232 -0.96133671  0.29628077 -0.20292438 -0.46606478
## [157] -0.82313394

6.4 Calculate ED Indexes

# Convert data to daily for P1_Luce
daily_date <- as.Date(time(precip_zoo))

# Daily precipitation totals for P1_Luce
daily_precip <- aggregate(precip_zoo$P1_Luce, daily_date, sum) 

# Check the basic statistics of the daily data
summary(daily_precip)

##      Index             daily_precip    
##  Min.   :2009-12-31   Min.   :  0.000  
##  1st Qu.:2013-03-31   1st Qu.:  0.000  
##  Median :2016-06-30   Median :  0.000  
##  Mean   :2016-06-30   Mean   :  4.236  
##  3rd Qu.:2019-09-30   3rd Qu.:  2.200  
##  Max.   :2022-12-30   Max.   :116.100

# Plot daily precipitation to examine data structure
plot(daily_precip, xlab="Year", ylab="Daily precipitation [mm]", main="Daily Precipitation for P1_Luce")

# Parameters
observation_days <- 365 
cumulative_precip <- rep(NA, length(daily_precip))

# Convert daily precipitation data to a numeric vector
precip_values <- as.numeric(daily_precip)  

# Calculate cumulative precipitation (EP)
for (day in (observation_days + 1):length(precip_values)) { 
  cumulative_sum <- 0
  for (past_day in 1:observation_days) {
    cumulative_sum <- cumulative_sum + sum(precip_values[(day - past_day + 1):day]) / past_day
  }
  cumulative_precip[day] <- cumulative_sum
}

# Mean of cumulative precipitation
mean_cumulative_precip <- mean(cumulative_precip, na.rm=TRUE)  

# Deviation from mean
deviation_precip <- cumulative_precip - mean_cumulative_precip 

# Standardized EDI
EDI <- deviation_precip / sd(deviation_precip, na.rm=TRUE)  

# Plot EDI with drought and wetness severity lines
plot(time(daily_precip), EDI, type="l", col="#000000", lty=3, 
     xlab="Year", ylab="EDI", main="Effective Drought Index (EDI) for P1_Luce", xaxt="n")

# Define yearly labels for the x-axis
years <- seq(2010, 2022, by = 1)
x_positions <- seq(from = as.numeric(time(daily_precip)[1]), 
                   to = as.numeric(time(daily_precip)[length(time(daily_precip))]), 
                   length.out = length(years))

# Add x-axis with yearly labels
axis(1, at = x_positions, labels = years)

# Add threshold lines
abline(h = c(-2, -1.5, -1, 1, 1.5, 2), 
       col = c("red", "orange", "yellow", "green", "lightblue", "blue"), 
       lty = 2)

# Add labels for each threshold
text(x = rep(mean(x_positions), 6), y = c(-2, -1.5, -1, 1, 1.5, 2), 
     labels = c("Extremely Dry", "Severely Dry", "Moderately Dry", 
                "Moderately Wet", "Severely Wet", "Extremely Wet"),
     col = c("red", "orange", "yellow", "green", "lightblue", "blue"), 
     pos = 4)

The plot shows the Effective Drought Index (EDI) values for P1_Luce, highlighting variations in drought and wetness conditions from 2010 to 2022. The years 2012 and mid-2021 display extremely dry conditions (EDI below -2), while mid-2012, late 2013, and mid-2014 indicate extremely wet conditions (EDI above 2). Seasonal fluctuations are evident, reflecting the influence of annual weather patterns.

# Plot SPI for comparison
par(mar = c(3, 3, 3, 3) + 0.1)
plot(time(monthly_precip), spi_3_P1_Luce$fitted, type="l", col="#000000", lty=3, 
     xlab="Year", ylab="SPI-3", main="Standardized Precipitation Index (SPI-3) for P1_Luce")

# Add threshold lines
abline(h = c(-2, -1.5, -1, 1, 1.5, 2), 
       col = c("red", "orange", "yellow", "green", "lightblue", "blue"), 
       lty = 2)

# Add labels for each threshold
text(x = rep(2011, 6), y = c(-2, -1.5, -1, 1, 1.5, 2),
     labels = c("Extremely Dry", "Severely Dry", "Moderately Dry", "Moderately Wet", "Severely Wet", "Extremely Wet"),
     col = c("red", "orange", "yellow", "green", "lightblue", "blue")) 

drought_zoo = zoo(selected$duration, selected$time)
par(new = TRUE)
plot(drought_zoo, axes = FALSE, xlab = "", ylab = "", type = "h", col = "purple", lwd = 2)

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "axes" is not a
## graphical parameter

axis(side = 4, at = pretty(range(drought_zoo)))
mtext("Duration of hydrological drought", side = 4, line = 3)

The plot compares the SPI-3 values for P1_Luce with the previously identified drought periods derived from the discharge data. Years such as 2010, 2011, 2012, 2017, and 2022 show consistency between meteorological droughts (negative SPI-3) and hydrological droughts (purple bars). However, some discrepancies are observed, such as in mid-2013, where precipitation data indicates wet conditions (positive SPI-3), while discharge data identifies drought conditions. This differences in meteorological and hydrological droughts are captured and emphasizes the complementary nature of these indices in understanding drought dynamics.

6.5 Correlation between NAO and precipitation data

# Download NAO index data
nao <- download_nao()

# Check the structure of the data
head(nao)  

## # A tibble: 6 × 3
##    Year Month   NAO
##   <int> <ord> <dbl>
## 1  1950 Jan    0.92
## 2  1950 Feb    0.4 
## 3  1950 Mar   -0.36
## 4  1950 Apr    0.73
## 5  1950 May   -0.59
## 6  1950 Jun   -0.06

# Basic statistics
summary(nao$NAO) 

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -3.180000 -0.750000  0.040000 -0.004378  0.720000  3.040000

# Plot NAO index data
plot(nao$NAO,ylab="NAO",type="l")

# Autocorrelation in NAO data
acf(nao$NAO[1:200])

# Aggregate Precipitation Data by Month
monthly_precip <- main_data %>%
  group_by(Year, Month) %>%
  summarize(
    P1_Luce = mean(P1_Luce, na.rm = TRUE),
    P2_Solcava = mean(P2_Solcava, na.rm = TRUE),
    P3_Lasko = mean(P3_Lasko, na.rm = TRUE),
    P4_GornjiGrad = mean(P4_GornjiGrad, na.rm = TRUE),
    P5_Radegunda = mean(P5_Radegunda, na.rm = TRUE)
  ) %>%
  mutate(YearMonth = as.yearmon(paste(Year, Month), "%Y %m"))

## `summarise()` has grouped output by 'Year'. You can override using the
## `.groups` argument.

# Prepare NAO Data
nao <- nao %>%
  mutate(YearMonth = as.yearmon(paste(Year, Month), "%Y %b"))

# Merge Precipitation and NAO Data
merged_data <- merge(monthly_precip, nao, by = "YearMonth", all.x = TRUE)

# Compute Correlations
correlations <- merged_data %>%
  dplyr::select(P1_Luce, P2_Solcava, P3_Lasko, P4_GornjiGrad, P5_Radegunda, NAO) %>%
  dplyr::summarize(across(-NAO, ~ cor(., merged_data$NAO))) %>%
  tidyr::pivot_longer(cols = everything(), names_to = "Station", values_to = "Correlation")


# Display correlations
correlations

## # A tibble: 5 × 2
##   Station       Correlation
##   <chr>               <dbl>
## 1 P1_Luce            -0.231
## 2 P2_Solcava         -0.228
## 3 P3_Lasko           -0.286
## 4 P4_GornjiGrad      -0.217
## 5 P5_Radegunda       -0.257

# Visualize correlations
ggplot(correlations, aes(x = Station, y = Correlation, fill = Station)) +
  geom_bar(stat = "identity") +
  labs(
    title = "Correlation between NAO and Precipitation",
    x = "Precipitation Station",
    y = "Correlation Coefficient"
  )

The bar chart show the correlation coefficients between the NAO (North Atlantic Oscillation) index and precipitation data across five meteorological stations. All stations exhibit a negative correlation, with values ranging from -0.217 at P4_GornjiGrad to -0.286 at P3_Lasko. This suggests that an increase in the NAO index is generally associated with decreased precipitation at these stations. The magnitude of these correlations highlights the varying sensitivity of different stations to atmospheric conditions influenced by the NAO.

7. SPATIAL DATA ANALYSIS

This section focuses on analyzing the spatial characteristics of the catchment. We will generate a stream network and catchment boundary using a digital terrain model and compare it with the actual stream network. Basic terrain statistics like slope will be calculated, along with a hypsometric curve to assess elevation distribution. Using the Thiessen polygon method, we will calculate areal precipitation and compare it to individual station data. Finally, land-use patterns will be analyzed using the 2012 Corine Land Cover (CLC) map to explore their impact on catchment hydrology.

7.1 Stream Network and Catchment Boundary Generation

# Import Catchment Area
catchment_area <- vect("basin/6210_Savinja_VSirje.shp")
catchment_area

##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 1, 2  (geometries, attributes)
##  extent      : 466905, 540775.8, 102666.7, 146214.3  (xmin, xmax, ymin, ymax)
##  source      : 6210_Savinja_VSirje.shp
##  coord. ref. : gauss_krueger_SLO 
##  names       : opomba      AREA
##  type        :  <chr>     <num>
##  values      :     NA 1.847e+09

# Import Digital Terrain Model (DTM)
digital_terrain_model <- rast("slo_100_ASCI.asc")
digital_terrain_model

## class       : SpatRaster 
## dimensions  : 1625, 2478, 1  (nrow, ncol, nlyr)
## resolution  : 100, 100  (x, y)
## extent      : 375212, 623012, 30872, 193372  (xmin, xmax, ymin, ymax)
## coord. ref. :  
## source      : slo_100_ASCI.asc 
## name        : slo_100_ASCI

# Plot Digital Terrain Model
plot(
  digital_terrain_model,
  main = "Digital Terrain Model (DEM)",
  col = terrain.colors(20),
  axes = TRUE
)

# Set the CRS for both the catchment area and the digital terrain model
crspost = crs(catchment_area)
crs(digital_terrain_model) = crspost

# Project the catchment area
project1 <- project(catchment_area, crspost)

# Plot the catchment area
plot(catchment_area, main = "Catchment Area Boundary", col = "lightblue", axes = TRUE)

# Overlay the projected catchment area
plot(project1, add = TRUE, border = "blue", lwd = 2)

# Add North Arrow and Scale Bar
north("topleft")
sbar(10000, xy = c(435000, 95000), type = "bar", divs = 2, below = "km", label = c(0, 5, 10))

7.2 Compare Generated Stream Network with Actual One

# Load Actual Stream Network
actual_stream_network <- vect("streams/VT_POVR_CRTLine.shp")


# Fill Local Depressions
wbt_breach_depressions_least_cost(dem = "slo_100_ASCI.asc", output = "dem_breached.tif", dist = 500, fill = TRUE)

# Fill Depressions Using Wang and Liu Algorithm
wbt_fill_depressions_wang_and_liu(dem = "dem_breached.tif", output = "filled_dem.tif")

# Flow Accumulation Calculation
wbt_d8_flow_accumulation(input = "filled_dem.tif", output = "flow_accumulation.tif")

# Flow Direction Calculation
wbt_d8_pointer(dem = "filled_dem.tif", output = "flow_direction.tif")

# Generate River Network
wbt_extract_streams(flow_accum = "flow_accumulation.tif", output = "generated_streams.tif", threshold = 1000)
                    

# Load Generated Stream Network
generated_stream_network <- rast("generated_streams.tif")

# Transform Generated Stream Network into Polygons
generated_stream_polygons <- as.polygons(generated_stream_network)

# Load Digital Terrain Model for Plotting
digital_terrain_model <- rast("slo_100_ASCI.asc")

# Compare Generated and Actual Stream Networks
plot(digital_terrain_model) 
lines(generated_stream_polygons, col = "gray")
lines(actual_stream_network, col = "red")

legend("topleft", legend = c("Generated Streams", "Actual Streams"), col = c("gray", "red"),
  lty = 1, lwd = 2)

plot(digital_terrain_model)

# Generate a catchment boundary
d1 = data.frame("veliko sirje", X=515395, Y=105435)
station = vect(d1, geom = c("X", "Y"))
points(station, lwd=100, col="red")

wbt_watershed(
  d8_pntr = "flow_direction.tif",
  pour_pts = "tocka.shp",
  output = "catchment_area.tif"
)

# Import the catchment area
catchment_area_raster <- rast("catchment_area.tif")

# Transform raster to polygon
catchment_area_polygon <- as.polygons(catchment_area_raster)

# Plot the digital terrain model
plot(digital_terrain_model, main = "Catchment Area and Digital Terrain Model")

# Add the catchment area to the plot
plot(catchment_area_polygon, add = TRUE, col = "orange", border = "red")
lines(catchment_area, col = "white")

7.3 Calculate Terrain Statistics for Selected Catchment

# Calculate Terrain Slope
terrain_slope <- terrain(digital_terrain_model, v = "slope", unit = "degrees")

# Calculate Basic Statistics
slope_statistics <- summary(values(terrain_slope))
slope_statistics

##      slope        
##  Min.   : 0.0     
##  1st Qu.: 3.4     
##  Median : 8.9     
##  Mean   :11.3     
##  3rd Qu.:16.7     
##  Max.   :72.0     
##  NA's   :1888211

# Visualize Slope
plot(terrain_slope, col = terrain.colors(10), main = "Slope Map")
north("topleft")

# Histogram of Slope
hist(values(terrain_slope), main = "Slope Histogram", xlab = "Slope (degrees)", breaks = 20,
     col = "lightblue")

The histogram represents the distribution of slope values within the Savinja Veliko Sirje catchment area. The majority of slopes are concentrated between 0° and 20°, with a significant peak in the lower range (0°–10°), indicating that the terrain is predominantly flat to gently sloping. Higher slope values (above 30°) are less frequent, suggesting that steep areas are limited in extent within the catchment. This slope distribution can have implications for runoff, erosion, and land-use planning within the region.

# Area and Perimeter of the Catchment
catchment_area_km2 <- expanse(catchment_area, unit = "km")
catchment_area_km2

## [1] 1845.736

catchment_perimeter <- perim(catchment_area)
catchment_perimeter

## [1] 272427.4

7.4 Plot Catchment Hypsometric Curve

# Plot Hypsometric Curve
plot(ecdf(as.numeric(values(digital_terrain_model))), xlab = "Elevation (m)",
     ylab = "Proportion of Area", main = "Hypsometric Curve of Catchment")

# Calculate and Display Elevation Quantiles
elevation_quantiles <- quantile(as.numeric(values(digital_terrain_model)), probs = c(0.25, 0.5, 0.75), na.rm = TRUE)
abline(v = elevation_quantiles, col = c("red", "blue", "green"), lty = 2, lwd = 2)

The hypsometric curve shows the distribution of elevation across the catchment area. The majority of the catchment’s area is concentrated at lower elevations, as indicated by the steep initial rise in the curve. Elevation quantiles are marked for reference: 25% of the catchment lies below approximately 288 meters (red line), 50% lies below 466 meters (blue line), and 75% lies below 719 meters (green line). This distribution suggests a predominance of lower elevation areas, which can influence runoff dynamics, sediment transport, and vegetation distribution in the catchment.

elevation_quantiles

##      25%      50%      75% 
## 287.6349 465.6508 719.4206

7.5 Thiessen Polygons and Areal Precipitation

# whitebox::install_whitebox()
precip_data <-  data.frame(names=c("P1_Luce", "P2_Solcava", "P3_Lasko", "P4_GornjiGrad", "P5_Radegunda"), lon=c(14.74, 14.69, 15.23, 14.8, 14.93), lat=c(46.35, 46.42, 46.15, 46.29, 46.36))

precip_data

##           names   lon   lat
## 1       P1_Luce 14.74 46.35
## 2    P2_Solcava 14.69 46.42
## 3      P3_Lasko 15.23 46.15
## 4 P4_GornjiGrad 14.80 46.29
## 5  P5_Radegunda 14.93 46.36

slo_model <- elevation_30s(country="SI", path="slo_100_ASCI")

xprecip <- vect(precip_data, geom=c("lon", "lat"))
crs(xprecip) <- crs(slo_model)
plot(slo_model)
points(xprecip, col = "red", lwd = 10)


theisson <- voronoi(xprecip)
lines(theisson, col = "green")

# reproject catchment boundary to the lat-lon geographic coordinate system
catchment_area <- vect("basin/6210_Savinja_VSirje.shp")
catchment_theisson <- project(catchment_area, crs(slo_model))
lines(catchment_theisson, col = "blue")


# Intersection between the catchment area and thiessen polygons
theisson_area <- intersect(catchment_theisson, theisson)
lines(theisson_area, lwd=2, col="red")

# calculate the areas of polygons 
expanse(theisson_area, unit = "km")

## [1] 160.6580 136.6530 923.7072 106.6205 518.0969

proportionP <- expanse(theisson_area, unit="km")/sum(expanse(theisson_area, unit="km"))
theisson_area$names

## [1] "P4_GornjiGrad" "P1_Luce"       "P3_Lasko"      "P2_Solcava"   
## [5] "P5_Radegunda"

# areal precipitation
areal_precip <- main_data$P1_Luce*proportionP[2] + main_data$P2_Solcava*proportionP[4] + main_data$P3_Lasko*proportionP[3] + main_data$P4_GornjiGrad*proportionP[1] + main_data$P5_Radegunda*proportionP[5]

main_data$Areal_Precip <- areal_precip
plot(main_data$Date, areal_precip, type = "l", xlab = "Date", ylab = "Areal Precipitation [mm]")

7.6 Analyse Land-Use of Catchment According to CLC Corine Map

# Load Land Cover map
land_use <- vect("land_use/CLC_2012_SIPolygon.shp")

# Plot the Land Cover map with the catchment area
plot(land_use, main = "Land Cover", axes = TRUE)
lines(catchment_area, col = "red")

# Check for invalid geometries in the land cover map and correct them
which(is.valid(land_use, TRUE) == FALSE)

##  [1] 11139 14587 14588 14613 14674 14680 14689 14720 14770 14816 14829 14912

land_use <- makeValid(land_use)

# Perform spatial intersection
intersations <- intersect(catchment_area, land_use)

head(values(intersations))

##   opomba       AREA        ID CODE_12   AREA_HA REMARK         OPIS3
## 1   <NA> 1847142640 SI-106579     231  48.46294   <NA>       Pašniki
## 2   <NA> 1847142640 SI-106609     231  48.02321   <NA>       Pašniki
## 3   <NA> 1847142640 SI-106655     311  88.05743   <NA> Listnati gozd
## 4   <NA> 1847142640 SI-106427     231 131.10541   <NA>       Pašniki
## 5   <NA> 1847142640 SI-106330     313 486.06507   <NA>   Mešani gozd
## 6   <NA> 1847142640 SI-106435     231  70.65305   <NA>       Pašniki
##                LABEL3   OPIS2   LABEL2
## 1            Pastures Pašniki Pastures
## 2            Pastures Pašniki Pastures
## 3 Broad-leaved forest Gozdovi  Forests
## 4            Pastures Pašniki Pastures
## 5        Mixed forest Gozdovi  Forests
## 6            Pastures Pašniki Pastures
##                                         OPIS1                        LABEL1
## 1                           Kmetijske povšine            Agricultural areas
## 2                           Kmetijske povšine            Agricultural areas
## 3 Gozdne in deloma ohranjene naravne površine Forest and semi natural areas
## 4                           Kmetijske povšine            Agricultural areas
## 5 Gozdne in deloma ohranjene naravne površine Forest and semi natural areas
## 6                           Kmetijske povšine            Agricultural areas
##   C_FAKTOR
## 1    0.004
## 2    0.004
## 3    0.002
## 4    0.004
## 5    0.002
## 6    0.004

# Plot the intersection
plot(intersations, "LABEL2", col = rainbow(20), plg = list(x = "bottomleft", cex = 0.7), mar = c(1, 1, 1, 1))

# Calculate the area of each land use
land_use_area = data.frame(area = expanse(intersations), desc = intersations$LABEL2)

# Calculate percentage for each land use type
land_use_area$percentage <- (land_use_area$area / sum(land_use_area$area)) * 100

head(land_use_area)

##          area     desc   percentage
## 1   54342.355 Pastures 0.0029442111
## 2   54005.764 Pastures 0.0029259750
## 3    7562.658  Forests 0.0004097368
## 4 1147967.363 Pastures 0.0621956535
## 5 3512183.833  Forests 0.1902863929
## 6  669435.041 Pastures 0.0362692801

# Group land use by land use type
land_use_summary <- land_use_area %>%
  group_by(desc) %>%
  summarise(
    total_area = sum(area),
    percentage = (sum(area) / sum(land_use_area$area)) * 100
  )

# Display summarized land use data
land_use_summary

## # A tibble: 12 × 3
##    desc                                             total_area percentage
##    <chr>                                                 <dbl>      <dbl>
##  1 Arable land                                       17402904.      0.943
##  2 Artificial, non-agricultural vegetated areas       2483004.      0.135
##  3 Forests                                         1032872691.     56.0  
##  4 Heterogeneous agricultural areas                 502325009.     27.2  
##  5 Industrial, commercial and transport units         9315112.      0.505
##  6 Inland waters                                      5256309.      0.285
##  7 Mine, dump and construction sites                  2329850.      0.126
##  8 Open spaces with little or no vegetation          20523663.      1.11 
##  9 Pastures                                         161496582.      8.75 
## 10 Permanent crops                                   17182187.      0.931
## 11 Scrub and/or herbaceous vegetation associations   39678600.      2.15 
## 12 Urban fabric                                      34869761.      1.89

land_use_colors <- rainbow(nrow(land_use_summary))

# Prepare data for plotting
land_use_summary$desc <- factor(land_use_summary$desc, levels = land_use_summary$desc)

# Create the plot
ggplot(land_use_summary, aes(x = percentage, y = desc, fill = desc)) +
  geom_bar(stat = "identity", color = "black", width = 0.5) + 
  scale_fill_manual(values = land_use_colors) + 
  labs(
    title = "Land Use Percentages in Catchment Area",
    x = "Percentage (%)",
    y = ""
  ) +
  theme_minimal() + 
  theme(
    legend.position = "none", 
    axis.text.y = element_text(size = 6, hjust = 1, color = "black"),  
    axis.text.x = element_text(size = 8),
    plot.title = element_text(size = 14, face = "bold")
  ) +
  geom_text(
    aes(label = paste0(round(percentage, 2), "%")), 
    hjust = -0.2, 
    color = "black",
    size = 3
  ) +
  coord_cartesian(clip = "off")

The bar chart represents the distribution of land use types within the catchment area. Forests dominate the catchment, covering 55.96% of the total area, followed by heterogeneous agricultural areas at 27.22%. Pastures account for 8.75%, while urban fabric and scrub or herbaceous vegetation associations represent 1.89% and 2.15%, respectively. Other land use types, including permanent crops, inland waters, and industrial or transport units, each occupy less than 1% of the catchment.

8. HYDROLOGICAL MODELLING (AIR GR)

In this section, we prepare, calibrate, validate, and compare hydrological models for the catchment. The GR4J and CemaNeige GR6J models will be used to simulate catchment hydrology. Relevant data will be prepared for both models, and datasets will be split (50% for calibration and validation) with the first year reserved as a warm-up period. Performance metrics of the models will be evaluated and compared to determine their suitability for modeling the catchment’s hydrology.

8.1 Prepare Data for GR4J and CemaNeige GR6J Models

# Create data frame
model_data <- main_data[, c("Date", "Discharge", "Runoff", "Areal_Precip", "Evapotrans", "Airtemp")]

# Display the new dataframe
head(model_data)

##         Date Discharge   Runoff Areal_Precip Evapotrans Airtemp
## 1 2010-01-01    59.058 2.762648    3.9654273        0.3     5.7
## 2 2010-01-02    53.925 2.522534    2.1476618        0.9     3.5
## 3 2010-01-03    47.681 2.230449    0.9421236        0.5    -2.0
## 4 2010-01-04    41.040 1.919792    0.0000000        0.2    -3.0
## 5 2010-01-05    38.694 1.810050    5.2556681        0.2    -1.8
## 6 2010-01-06    36.966 1.729216    7.3304657        0.2    -1.0

# Display basic statistics
summary(model_data)

##       Date                          Discharge           Runoff       
##  Min.   :2010-01-01 00:00:00.00   Min.   :  5.515   Min.   : 0.2580  
##  1st Qu.:2013-04-01 18:00:00.00   1st Qu.: 13.966   1st Qu.: 0.6533  
##  Median :2016-07-01 12:00:00.00   Median : 22.032   Median : 1.0306  
##  Mean   :2016-07-01 12:24:47.61   Mean   : 38.480   Mean   : 1.8000  
##  3rd Qu.:2019-10-01 06:00:00.00   3rd Qu.: 40.252   3rd Qu.: 1.8829  
##  Max.   :2022-12-31 00:00:00.00   Max.   :859.744   Max.   :40.2176  
##   Areal_Precip         Evapotrans       Airtemp      
##  Min.   :  0.00000   Min.   :0.000   Min.   :-13.80  
##  1st Qu.:  0.00000   1st Qu.:0.700   1st Qu.:  4.40  
##  Median :  0.01672   Median :1.800   Median : 11.30  
##  Mean   :  3.60893   Mean   :2.264   Mean   : 10.88  
##  3rd Qu.:  2.57378   3rd Qu.:3.600   3rd Qu.: 17.70  
##  Max.   :109.52009   Max.   :7.700   Max.   : 28.90

# Plot flow data
plot(model_data$Runoff, type = "l", col = "blue", lwd = 2, main = "Flow Data Over Time", xlab = "Days", ylab = "Flow (mm)")

8.2 Calibrate and Validate Both Hydrological Models

# Create Inputs Model
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = as.POSIXct(model_data$Date, format = "%Y-%m-%d"), Precip = model_data$Areal_Precip, PotEvap = model_data$Evapotrans)

# Identify the index where calibration ends
calibration_end <- which(format(as.Date(model_data$Date), format = "%Y-%m-%d") == "2016-12-31")

# Define the sequence for the run period
Ind_Run <- seq(366, calibration_end)


# Create RunOptions with defined run and warm-up periods
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel, IndPeriod_Run = Ind_Run, IndPeriod_WarmUp = 1:365)

# Extract observed discharge data for the run period
Obs_Discharge <- model_data$Runoff[Ind_Run]

# Create InputsCrit using the Nash-Sutcliffe Efficiency (NSE) criterion
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel, RunOptions = RunOptions, Obs = Obs_Discharge)

# Define calibration options using the Michel algorithm
CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibration_Michel)

# Perform calibration to find optimal model parameters
OutputsCalib <- Calibration_Michel(InputsModel = InputsModel, RunOptions = RunOptions, InputsCrit = InputsCrit, CalibOptions = CalibOptions, FUN_MOD = RunModel_GR4J)

## Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
##   Screening completed (81 runs)
##       Param =  432.681,   -0.020,   42.098,    1.417
##       Crit. NSE[Q]       = 0.5850
## Steepest-descent local search in progress
##   Calibration completed (27 iterations, 282 runs)
##       Param =  549.741,    0.270,   45.955,    0.548
##       Crit. NSE[Q]       = 0.8192

# Extract calibrated parameter values
ParamGR4J <- OutputsCalib$ParamFinalR

ParamGR4J 

## [1] 549.7413151   0.2699136  45.9554779   0.5477049

# Run the GR4J model with calibrated parameters
OutputsModel <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions, Param = ParamGR4J)

# Plot the comparison between measured and modeled data
plot(OutputsModel, model_data$Runoff[366:calibration_end])

# Define the sequence for the validation period
Ind_Run1 <- seq((calibration_end + 1), nrow(model_data))

# Create RunOptions for the validation period
RunOptions1 <- CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel, IndPeriod_Run = Ind_Run1)

## Warning in CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel, : model warm up period not defined: default configuration used
##   the year preceding the run period is used

# Run the GR4J model for the validation period
OutputsModel1 <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions1, Param = ParamGR4J)

# Plot the comparison between measured and modeled runoff for validation
plot(OutputsModel1, model_data$Runoff[Ind_Run1])

# Create InputsCrit for validation using NSE as the criterion
InputsCrit1 <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel, RunOptions = RunOptions1, Obs = model_data$Runoff[Ind_Run1])

# Calculate the validation performance (NSE)
OutputsCrit1 <- ErrorCrit_NSE(InputsCrit = InputsCrit1, OutputsModel = OutputsModel1)

## Crit. NSE[Q] = 0.7192

# Print the validation results
OutputsCrit1

## $CritValue
## [1] 0.7191667
## 
## $CritName
## [1] "NSE[Q]"
## 
## $CritBestValue
## [1] 1
## 
## $Multiplier
## [1] -1
## 
## $Ind_notcomputed
## NULL
## 
## attr(,"class")
## [1] "NSE"       "ErrorCrit"

# Prepare observation data in the required format
preob <- data.frame(DatesR = as.POSIXct(model_data$Date, format = "%Y-%m-%d"), P=model_data[,4], E=model_data[, 5],Qmm=model_data[,3]) 

# Display the first few rows of the observation data
head(preob)

##       DatesR         P   E      Qmm
## 1 2010-01-01 3.9654273 0.3 2.762648
## 2 2010-01-02 2.1476618 0.9 2.522534
## 3 2010-01-03 0.9421236 0.5 2.230449
## 4 2010-01-04 0.0000000 0.2 1.919792
## 5 2010-01-05 5.2556681 0.2 1.810050
## 6 2010-01-06 7.3304657 0.2 1.729216

# Load example data for the basin
data(L0123001) 

# Display the hypsometric curve data format
BasinInfo

## $BasinCode
## [1] "L0123001"
## 
## $BasinName
## [1] "Blue River at Nourlangie Rock"
## 
## $BasinArea
## [1] 360
## 
## $HypsoData
##   [1]  286  309  320  327  333  338  342  347  351  356  360  365  369  373  378
##  [16]  382  387  393  399  405  411  417  423  428  434  439  443  448  453  458
##  [31]  463  469  474  480  485  491  496  501  507  513  519  524  530  536  542
##  [46]  548  554  560  566  571  577  583  590  596  603  609  615  622  629  636
##  [61]  642  649  656  663  669  677  684  691  698  706  714  722  730  738  746
##  [76]  754  762  770  777  786  797  808  819  829  841  852  863  875  887  901
##  [91]  916  934  952  972  994 1012 1029 1054 1080 1125 1278

# Calculate the hypsometric curve using the digital terrain model
hypso <- c(min(as.numeric(values(digital_terrain_model)),na.rm=T),as.numeric(quantile(as.numeric(values(digital_terrain_model)),probs=seq(from=0.01,by=0.01,to=0.99),na.rm=T)),max(as.numeric(values(digital_terrain_model)),na.rm=T))

# Define model settings, including the number of height zones and the hypsometric curve
InputsModel2 <- CreateInputsModel(FUN_MOD = RunModel_CemaNeigeGR6J, as.POSIXct(model_data$Date, format = "%Y-%m-%d"), Precip = model_data[,4], PotEvap = model_data[, 5], TempMean = model_data[,6],  ZInputs = median(hypso), HypsoData = hypso, NLayers = 5)

## input series were successfully created on 5 elevation layers for use by CemaNeige

# Define run options, including warm-up period and annual solid precipitation
RunOptions2 <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR6J, InputsModel = InputsModel2, IndPeriod_Run = Ind_Run, IndPeriod_WarmUp = 1:365, MeanAnSolidPrecip=c(70,130,180,250,350)) 

# Define the goodness of fit criterion and discharge data
InputsCrit2 <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel2, RunOptions = RunOptions2, 
                                Obs = model_data[366:calibration_end,3]) 

# Define the calibration method
CalibOptions2 <- CreateCalibOptions(FUN_MOD = RunModel_CemaNeigeGR6J, FUN_CALIB = Calibration_Michel) 

# Run the model calibration
OutputsCalib2 <- Calibration_Michel(InputsModel = InputsModel2, RunOptions = RunOptions2, InputsCrit = InputsCrit2, CalibOptions = CalibOptions2, FUN_MOD = RunModel_CemaNeigeGR6J) 

## Grid-Screening in progress (

## 0% 20% 40% 60% 80% 100%)
##   Screening completed (6561 runs)
##       Param =   90.017,   -0.521,   27.113,    1.369,    0.220,   54.598,    0.002,    6.764
##       Crit. NSE[Q]       = 0.6486
## Steepest-descent local search in progress
##   Calibration completed (65 iterations, 7565 runs)
##       Param =  178.520,    0.107,   14.540,    1.042,    0.447,   32.407,    0.002,    6.540
##       Crit. NSE[Q]       = 0.8876

# Save the calibrated parameter values
ParamGR6JCemaNeige <- OutputsCalib2$ParamFinalR 

# Run the model with the calibrated parameter values
OutputsModel2 <- RunModel_CemaNeigeGR6J(InputsModel = InputsModel2, RunOptions = RunOptions2, 
                                        Param = ParamGR6JCemaNeige) 

# Plot the observed versus simulated values
plot(OutputsModel2, model_data[366:calibration_end,3]) 

# Calculate NSE for the validation period
ErrorCrit_NSE(InputsCrit = InputsCrit2, OutputsModel = OutputsModel2) 

## Crit. NSE[Q] = 0.8876

## $CritValue
## [1] 0.8876181
## 
## $CritName
## [1] "NSE[Q]"
## 
## $CritBestValue
## [1] 1
## 
## $Multiplier
## [1] -1
## 
## $Ind_notcomputed
## NULL
## 
## attr(,"class")
## [1] "NSE"       "ErrorCrit"

# Create run options for validation
RunOptions3 <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR6J, InputsModel = InputsModel2, IndPeriod_Run = Ind_Run1) 

## Warning in CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR6J, InputsModel = InputsModel2, : model warm up period not defined: default configuration used
##   the year preceding the run period is used

## Warning in CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR6J, InputsModel = InputsModel2, : 'MeanAnSolidPrecip' not defined: it was automatically set to c(146,146,146,146,146)  from the 'InputsModel' given to the function. Be careful in case your application is (short-term) forecasting.

# Run the model for validation
OutputsModel3 <- RunModel_CemaNeigeGR6J(InputsModel = InputsModel2, RunOptions = RunOptions3, Param = ParamGR6JCemaNeige)

# Plot simulated vs observed discharge
plot(OutputsModel3, Qobs = model_data[Ind_Run1,3]) 

# Create inputs for the NSE criterion during validation
InputsCrit3 <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel2, RunOptions = RunOptions3, Obs = model_data[Ind_Run1,3]) 

# Calculate NSE for the validation period
OutputsCrit3 <- ErrorCrit_NSE(InputsCrit = InputsCrit3, OutputsModel = OutputsModel3) 

## Crit. NSE[Q] = 0.7892

8.3 Compare the Performance of GR4J and GR6J Models

# Calculate NSE for GR4J model
OutputsCrit1 <- ErrorCrit_NSE(InputsCrit = InputsCrit1, OutputsModel = OutputsModel1)

## Crit. NSE[Q] = 0.7192

NSE_GR4J <- OutputsCrit1$CritValue 

# Calculate NSE for GR6J model
OutputsCrit3 <- ErrorCrit_NSE(InputsCrit = InputsCrit3, OutputsModel = OutputsModel3)

## Crit. NSE[Q] = 0.7892

NSE_GR6J <- OutputsCrit3$CritValue 

# Compare the performance of the two models
if (NSE_GR6J > NSE_GR4J) {
  cat("The GR6J model performs better with NSE =", NSE_GR6J,
      "compared to GR4J with NSE =", NSE_GR4J, "\n")
} else if (NSE_GR6J < NSE_GR4J) {
  cat("The GR4J model performs better with NSE =", NSE_GR4J,
      "compared to GR6J with NSE =", NSE_GR6J, "\n")
} else {
  cat("Both models perform equally well with NSE =", NSE_GR6J, "\n")
}

## The GR6J model performs better with NSE = 0.7891693 compared to GR4J with NSE = 0.7191667

The comparison of model performance based on Nash-Sutcliffe Efficiency (NSE) indicates that the GR6J model outperforms the GR4J model. The GR6J model achieved an NSE value of 0.789, whereas the GR4J model achieved an NSE value of 0.719. The enhanced performance of the GR6J model can be attributed to its additional parameters, including an exponential store that improves low-flow simulations, allowing it to better capture hydrological processes compared to the GR4J model. This suggests that GR6J provides a more accurate representation of observed streamflow dynamics.

9. HYDROLOGICAL MODELLING (TUWIEN MODEL)

In this section, we will explore hydrological modeling using the TUWIEN model. The analysis begins by preparing all relevant data for the model. Following data preparation, the TUWIEN model will be calibrated and validated by splitting the dataset into training and validation periods, using the first year of data for warm-up. Finally, the performance of the TUWIEN model will be compared with the CemaNeige GR6J model. Both graphical and statistical metrics will be used to assess the models’ capabilities.

9.1 Prepare data for the TUWIEN model

# Create a data frame for the TUWIEN model
preob1 <- data.frame(
  DatesR = as.POSIXct(model_data$Date, format = "%Y-%m-%d"),
  P = model_data[, 4],
  E = model_data[, 5],
  Qmm = model_data[, 3],
  T = model_data[, 6]
)
head(preob1)

##       DatesR         P   E      Qmm    T
## 1 2010-01-01 3.9654273 0.3 2.762648  5.7
## 2 2010-01-02 2.1476618 0.9 2.522534  3.5
## 3 2010-01-03 0.9421236 0.5 2.230449 -2.0
## 4 2010-01-04 0.0000000 0.2 1.919792 -3.0
## 5 2010-01-05 5.2556681 0.2 1.810050 -1.8
## 6 2010-01-06 7.3304657 0.2 1.729216 -1.0

# Plot precipitation and runoff data
plot(preob1$DatesR, preob1$P, type = "l", col = "blue", ylab = "Precipitation [mm]", xlab = "Date", main = "Precipitation over Time")

plot(preob1$DatesR, preob1$Qmm, type = "l", col = "green", ylab = "Runoff [mm]", xlab = "Date", main = "Runoff over Time")

9.2 Calibrate and validate TUWIEN model

# End of calibration period
konec <- which(format(preob1$DatesR, format = "%Y-%m-%d") == "2016-12-31")

# Run TUWIEN model with randomly selected parameter values
sim1 <- TUWmodel(
  prec = preob1$P[1:konec],
  airt = preob1$T[1:konec],
  ep = preob1$E[1:konec],
  area = 1,
  param = c(1.3, 2.0, -1.0, 1.0, 0.0, 0.8, 360.0, 0.2, 0.3, 7.0, 150.0, 50.0, 2.0, 10.0, 25.0),
  incon = c(50, 0, 2.5, 2.5)
)

# Check the structure of the simulation output
str(sim1)

## List of 22
##  $ itsteps: int 2557
##  $ nzones : int 1
##  $ area   : num 1
##  $ param  : num [1, 1:15] 1.3 2 -1 1 0 0.8 360 0.2 0.3 7 ...
##   ..- attr(*, "names")= chr [1:15] "SCF" "DDF" "Tr" "Ts" ...
##  $ incon  : num [1, 1:4] 50 0 2.5 2.5
##   ..- attr(*, "names")= chr [1:4] "SSM0" "SWE0" "SUZ0" "SLZ0"
##  $ prec   : num [1:2557] 3.965 2.148 0.942 0 5.256 ...
##  $ airt   : num [1:2557] 5.7 3.5 -2 -3 -1.8 -1 -0.3 -0.4 0.6 0.9 ...
##  $ ep     : num [1:2557] 0.3 0.9 0.5 0.2 0.2 0.2 0.2 0.2 0.2 0.2 ...
##  $ output : num [1, 1:20, 1:2557] 0.404 0 51.24 3.965 0 ...
##  $ qzones : num [1, 1:2557] 0.40409 0.1516 0.15334 0.00696 0.01563 ...
##  $ q      : num [1, 1:2557] 0.40409 0.1516 0.15334 0.00696 0.01563 ...
##  $ swe    : num [1, 1:2557] 0 0 1.22 1.22 8.06 ...
##  $ melt   : num [1, 1:2557] 0 0 0 0 0 0 0 0 1.2 1.8 ...
##  $ q0     : num [1, 1:2557] 0 0 0 0 0 0 0 0 0 0 ...
##  $ q1     : num [1, 1:2557] 0.374 0.26 0 0 0 ...
##  $ q2     : num [1, 1:2557] 0.0298 0.0429 0.0558 0.0554 0.055 ...
##  $ moist  : num [1, 1:2557] 51.2 51.8 51.8 51.8 51.8 ...
##  $ rain   : num [1, 1:2557] 3.97 2.15 0 0 0 ...
##  $ snow   : num [1, 1:2557] 0 0 0.942 0 5.256 ...
##  $ eta    : num [1, 1:2557] 0.0534 0.1623 0 0 0 ...
##  $ suz    : num [1, 1:2557] 2.8 1.99 0 0 0 ...
##  $ slz    : num [1, 1:2557] 4.47 6.43 8.36 8.31 8.25 ...

# Plot simulated and observed runoff for the calibration period
plot(as.numeric(sim1$q)[366:konec], type = "l", col = "red", ylab = "Discharge [mm]", xlab = "Days", main = "Calibration Results")
lines(preob1$Qmm[366:konec], col = "green")
legend("topright", legend = c("Simulated", "Observed"), col = c("red", "green"), lty = 1)

# plot the simulated values (for the first year)
plot(as.numeric(sim1$q)[1:365],type="l", col="red",ylab="Discharge [mm]") 

# add measured flow data
lines(preob1[1:365,4], col="green") 

# Calculate goodness-of-fit statistics for the calibration period
gof_results <- gof(as.numeric(sim1$q)[366:konec], preob1$Qmm[366:konec])
gof_results

##          [,1]
## ME       1.11
## MAE      1.45
## MSE      3.67
## RMSE     1.92
## ubRMSE   1.56
## NRMSE % 76.30
## PBIAS % 59.80
## RSR      0.76
## rSD      0.86
## NSE      0.42
## mNSE    -0.02
## rNSE    -0.74
## wNSE     0.74
## wsNSE    0.27
## d        0.83
## dr       0.49
## md       0.50
## rd       0.49
## cp      -0.09
## r        0.79
## R2       0.42
## bR2      0.42
## VE       0.22
## KGE      0.35
## KGElf    0.20
## KGEnp    0.33
## KGEkm    0.32

# Define the objective function for optimization
msespr <- function (param, precip, temp, potevap, runoff, area) { 
  simu <- TUWmodel(param, prec=as.numeric(precip), airt=as.numeric(temp), ep=as.numeric(potevap), area=area)$q 
  simu <- simu[-c(1:365)]  # Remove warm-up period
  obse <- runoff[-c(1:365)]
  mse(simu, obse)  # Calculate mean squared error
}

# Calibrate the model using DEoptim
calibrate_period1 <- DEoptim(
  fn = msespr,
  lower = c(0.9, 0.0, 1.0, -3.0, -2.0, 0.0, 0.0, 0.0, 0.0, 2.0, 30.0, 1.0, 0.0, 0.0, 0.0),
  upper = c(1.5, 5.0, 3.0, 1.0, 2.0, 1.0, 600.0, 20.0, 2.0, 30.0, 250.0, 100.0, 8.0, 30.0, 50.0),
  control = DEoptim.control(itermax = 60),
  precip = preob1$P[1:konec],
  temp = preob1$T[1:konec],
  potevap = preob1$E[1:konec],
  runoff = preob1$Qmm[1:konec],
  area = 1
)

## Iteration: 1 bestvalit: 1.562010 bestmemit:    1.258932    0.834634    2.819617   -0.828915   -1.962681    0.986918  319.547874   17.270585    1.128826    7.168388  177.013609   22.856899    3.719463    2.875780   40.803773
## Iteration: 2 bestvalit: 1.562010 bestmemit:    1.258932    0.834634    2.819617   -0.828915   -1.962681    0.986918  319.547874   17.270585    1.128826    7.168388  177.013609   22.856899    3.719463    9.268427   40.803773
## Iteration: 3 bestvalit: 1.515548 bestmemit:    1.134806    4.458489    2.786875   -2.867137    0.428393    0.875656  240.440627   19.381912    0.676021   26.449547  107.558175    3.290968    6.649246   23.730181   17.794380
## Iteration: 4 bestvalit: 1.453108 bestmemit:    1.053977    2.367076    1.982698   -1.581068   -1.023994    0.672748  105.929010   17.989298    1.721246    7.707436  168.114704   29.029833    3.775845   16.283480   11.133193
## Iteration: 5 bestvalit: 1.264874 bestmemit:    1.105428    0.707724    1.829851    0.186863   -0.934620    0.946221  279.553487    8.838589    1.762308    8.537585   38.952472    4.706629    5.459383   23.533203   46.675108
## Iteration: 6 bestvalit: 1.264874 bestmemit:    1.105428    0.707724    1.829851    0.186863   -0.934620    0.946221  279.553487    8.838589    1.762308    8.537585   38.952472    4.706629    5.459383   23.533203   46.675108
## Iteration: 7 bestvalit: 1.236245 bestmemit:    1.430549    0.786743    2.201725   -2.233606   -0.381786    0.998964  352.814081    2.247926    1.813976    2.156076   40.920631   14.606405    2.894247   10.711176   49.275859
## Iteration: 8 bestvalit: 1.231692 bestmemit:    1.430549    0.786743    2.006815   -2.233606   -0.381786    0.998964  352.814081    2.247926    1.813976    2.156076   40.920631   14.606405    2.894247   10.711176   49.275859
## Iteration: 9 bestvalit: 1.231692 bestmemit:    1.430549    0.786743    2.006815   -2.233606   -0.381786    0.998964  352.814081    2.247926    1.813976    2.156076   40.920631   14.606405    2.894247   10.711176   49.275859
## Iteration: 10 bestvalit: 1.198061 bestmemit:    0.943814    1.478951    1.403058   -1.183636   -1.956553    0.700740  315.402509    1.717766    1.238394    3.304905  151.668674   23.085526    1.803697   28.368848   49.664354
## Iteration: 11 bestvalit: 1.198061 bestmemit:    0.943814    1.478951    1.403058   -1.183636   -1.956553    0.700740  315.402509    1.717766    1.238394    3.304905  151.668674   23.085526    1.803697   28.368848   49.664354
## Iteration: 12 bestvalit: 1.198061 bestmemit:    0.943814    1.478951    1.403058   -1.183636   -1.956553    0.700740  315.402509    1.717766    1.238394    3.304905  151.668674   23.085526    1.803697   28.368848   49.664354
## Iteration: 13 bestvalit: 1.198061 bestmemit:    0.943814    1.478951    1.403058   -1.183636   -1.956553    0.700740  315.402509    1.717766    1.238394    3.304905  151.668674   23.085526    1.803697   28.368848   49.664354
## Iteration: 14 bestvalit: 1.198061 bestmemit:    0.943814    1.478951    1.403058   -1.183636   -1.956553    0.700740  315.402509    1.717766    1.238394    3.304905  151.668674   23.085526    1.803697   28.368848   49.664354
## Iteration: 15 bestvalit: 0.985166 bestmemit:    0.948157    1.080407    2.006815   -2.233606   -0.381786    0.998964  352.814081    2.247926    1.813976    2.156076   40.920631   14.606405    2.894247   10.711176   49.275859
## Iteration: 16 bestvalit: 0.985166 bestmemit:    0.948157    1.080407    2.006815   -2.233606   -0.381786    0.998964  352.814081    2.247926    1.813976    2.156076   40.920631   14.606405    2.894247    0.479086   32.196982
## Iteration: 17 bestvalit: 0.985166 bestmemit:    0.948157    1.080407    2.006815   -2.233606   -0.381786    0.998964  352.814081    2.247926    1.813976    2.156076   40.920631   14.606405    2.894247    0.479086   32.196982
## Iteration: 18 bestvalit: 0.985061 bestmemit:    0.948157    1.080407    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.813976    2.156076   40.920631   14.606405    2.894247    0.479086   32.196982
## Iteration: 19 bestvalit: 0.985061 bestmemit:    0.948157    1.080407    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.813976    2.156076   40.920631   14.606405    2.894247    0.479086   32.196982
## Iteration: 20 bestvalit: 0.985061 bestmemit:    0.948157    1.080407    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.813976    2.156076   40.920631   14.606405    2.894247    0.479086   32.196982
## Iteration: 21 bestvalit: 0.985061 bestmemit:    0.948157    1.080407    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.813976    2.156076   40.920631   14.606405    2.894247    0.479086   32.196982
## Iteration: 22 bestvalit: 0.985061 bestmemit:    0.948157    1.080407    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.813976    2.156076   40.920631   14.606405    2.894247    0.479086   32.196982
## Iteration: 23 bestvalit: 0.965310 bestmemit:    0.948157    1.080407    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.813976    2.156076   40.920631   16.587395    2.894247    0.479086   32.196982
## Iteration: 24 bestvalit: 0.949514 bestmemit:    0.948157    1.392848    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.813976    2.156076   40.920631   16.587395    2.894247    0.479086   32.196982
## Iteration: 25 bestvalit: 0.949514 bestmemit:    0.948157    1.392848    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.813976    2.156076   40.920631   16.587395    2.894247    0.479086   32.196982
## Iteration: 26 bestvalit: 0.949514 bestmemit:    0.948157    1.392848    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.813976    2.156076   40.920631   16.587395    2.894247    0.479086   32.196982
## Iteration: 27 bestvalit: 0.949514 bestmemit:    0.948157    1.392848    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.813976    2.156076   40.920631   16.587395    2.894247    0.479086   32.196982
## Iteration: 28 bestvalit: 0.949514 bestmemit:    0.948157    1.392848    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.813976    2.156076   40.920631   16.587395    2.894247    0.479086   32.196982
## Iteration: 29 bestvalit: 0.949514 bestmemit:    0.948157    1.392848    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.813976    2.156076   40.920631   16.587395    2.894247    0.479086   32.196982
## Iteration: 30 bestvalit: 0.949514 bestmemit:    0.948157    1.392848    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.813976    2.156076   40.920631   16.587395    2.894247    0.479086   32.196982
## Iteration: 31 bestvalit: 0.949514 bestmemit:    0.948157    1.392848    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.813976    2.156076   40.920631   16.587395    2.894247    0.479086   32.196982
## Iteration: 32 bestvalit: 0.949514 bestmemit:    0.948157    1.392848    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.813976    2.156076   40.920631   16.587395    2.894247    0.479086   32.196982
## Iteration: 33 bestvalit: 0.949514 bestmemit:    0.948157    1.392848    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.813976    2.156076   40.920631   16.587395    2.894247    0.479086   32.196982
## Iteration: 34 bestvalit: 0.946884 bestmemit:    0.948157    1.392848    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.932069    2.156076   40.920631   16.587395    2.894247    0.479086   32.196982
## Iteration: 35 bestvalit: 0.946884 bestmemit:    0.948157    1.392848    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.932069    2.156076   40.920631   16.587395    2.894247    0.479086   32.196982
## Iteration: 36 bestvalit: 0.946884 bestmemit:    0.948157    1.392848    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.932069    2.156076   40.920631   16.587395    2.894247    0.479086   32.196982
## Iteration: 37 bestvalit: 0.946884 bestmemit:    0.948157    1.392848    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.932069    2.156076   40.920631   16.587395    2.894247    0.479086   32.196982
## Iteration: 38 bestvalit: 0.946884 bestmemit:    0.948157    1.392848    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.932069    2.156076   40.920631   16.587395    2.894247    0.479086   32.196982
## Iteration: 39 bestvalit: 0.946884 bestmemit:    0.948157    1.392848    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.932069    2.156076   40.920631   16.587395    2.894247    0.479086   32.196982
## Iteration: 40 bestvalit: 0.946884 bestmemit:    0.948157    1.392848    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.932069    2.156076   40.920631   16.587395    2.894247    0.479086   32.196982
## Iteration: 41 bestvalit: 0.946884 bestmemit:    0.948157    1.392848    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.932069    2.156076   40.920631   16.587395    2.894247    0.479086   32.196982
## Iteration: 42 bestvalit: 0.936877 bestmemit:    0.922315    1.392848    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.932069    2.156076   40.920631   16.587395    2.894247    0.479086   32.196982
## Iteration: 43 bestvalit: 0.936877 bestmemit:    0.922315    1.392848    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.932069    2.156076   40.920631   16.587395    2.894247    0.479086   32.196982
## Iteration: 44 bestvalit: 0.936877 bestmemit:    0.922315    1.392848    2.006815   -2.233606   -0.269818    0.998964  352.814081    2.247926    1.932069    2.156076   40.920631   16.587395    2.894247    0.479086   32.196982
## Iteration: 45 bestvalit: 0.935961 bestmemit:    0.960670    1.386594    2.115617    0.171584   -1.342859    0.958450  297.556700    3.363709    1.763296    7.986158   34.336170    6.521182    4.303730    8.403318   33.161193
## Iteration: 46 bestvalit: 0.935961 bestmemit:    0.960670    1.386594    2.115617    0.171584   -1.342859    0.958450  297.556700    3.363709    1.763296    7.986158   34.336170    6.521182    4.303730    8.403318   33.161193
## Iteration: 47 bestvalit: 0.935961 bestmemit:    0.960670    1.386594    2.115617    0.171584   -1.342859    0.958450  297.556700    3.363709    1.763296    7.986158   34.336170    6.521182    4.303730    8.403318   33.161193
## Iteration: 48 bestvalit: 0.935961 bestmemit:    0.960670    1.386594    2.115617    0.171584   -1.342859    0.958450  297.556700    3.363709    1.763296    7.986158   34.336170    6.521182    4.303730    8.403318   33.161193
## Iteration: 49 bestvalit: 0.935961 bestmemit:    0.960670    1.386594    2.115617    0.171584   -1.342859    0.958450  297.556700    3.363709    1.763296    7.986158   34.336170    6.521182    4.303730    8.403318   33.161193
## Iteration: 50 bestvalit: 0.935961 bestmemit:    0.960670    1.386594    2.115617    0.171584   -1.342859    0.958450  297.556700    3.363709    1.763296    7.986158   34.336170    6.521182    4.303730    8.403318   33.161193
## Iteration: 51 bestvalit: 0.935961 bestmemit:    0.960670    1.386594    2.115617    0.171584   -1.342859    0.958450  297.556700    3.363709    1.763296    7.986158   34.336170    6.521182    4.303730    8.403318   33.161193
## Iteration: 52 bestvalit: 0.935961 bestmemit:    0.960670    1.386594    2.115617    0.171584   -1.342859    0.958450  297.556700    3.363709    1.763296    7.986158   34.336170    6.521182    4.303730    8.403318   33.161193
## Iteration: 53 bestvalit: 0.935961 bestmemit:    0.960670    1.386594    2.115617    0.171584   -1.342859    0.958450  297.556700    3.363709    1.763296    7.986158   34.336170    6.521182    4.303730    8.403318   33.161193
## Iteration: 54 bestvalit: 0.935121 bestmemit:    0.960670    1.386594    1.892064    0.171584   -1.342859    0.958450  297.556700    3.363709    1.763296    7.986158   34.336170    6.521182    4.303730    8.403318   33.161193
## Iteration: 55 bestvalit: 0.933865 bestmemit:    0.960670    1.386594    1.892064    0.330761   -1.342859    0.958450  297.556700    3.363709    1.763296    7.986158   34.336170    6.521182    4.303730    8.403318   33.161193
## Iteration: 56 bestvalit: 0.933865 bestmemit:    0.960670    1.386594    1.892064    0.330761   -1.342859    0.958450  297.556700    3.363709    1.763296    7.986158   34.336170    6.521182    4.303730    8.403318   33.161193
## Iteration: 57 bestvalit: 0.895191 bestmemit:    0.960670    1.386594    1.892064    0.330761   -1.342859    0.958450  297.556700    3.363709    1.763296    7.986158   34.336170    3.314976    4.303730    8.403318   33.161193
## Iteration: 58 bestvalit: 0.895191 bestmemit:    0.960670    1.386594    1.892064    0.330761   -1.342859    0.958450  297.556700    3.363709    1.763296    7.986158   34.336170    3.314976    4.303730    8.403318   33.161193
## Iteration: 59 bestvalit: 0.895191 bestmemit:    0.960670    1.386594    1.892064    0.330761   -1.342859    0.958450  297.556700    3.363709    1.763296    7.986158   34.336170    3.314976    4.303730    8.403318   33.161193
## Iteration: 60 bestvalit: 0.895191 bestmemit:    0.960670    1.386594    1.892064    0.330761   -1.342859    0.958450  297.556700    3.363709    1.763296    7.986158   34.336170    3.314976    4.303730    8.403318   33.161193

# Run the model with calibrated parameters
sim2 <- TUWmodel(
  prec = preob1$P[1:konec],
  airt = preob1$T[1:konec],
  ep = preob1$E[1:konec],
  area = 1,
  param = calibrate_period1$optim$bestmem
)

# Plot simulated vs observed runoff for the validation period
plot(as.numeric(sim2$q)[366:konec], type = "l", col = "red", ylab = "Discharge [mm]", xlab = "Days", main = "Validation Results")
lines(preob1$Qmm[366:konec], col = "green")
legend("topright", legend = c("Simulated", "Observed"), col = c("red", "green"), lty = 1)

# Plot simulated vs observed runoff for the first year
plot(as.numeric(sim2$q)[1:365], type = "l", col = "red", ylab = "Discharge [mm]", xlab = "Days", main = "Validation Results for First Year") 
lines(preob1[1:365, 4], col = "green")

# Calculate NSE for the validation period
NSE_validation <- NSE(as.numeric(sim2$q)[366:konec], preob1$Qmm[366:konec])
print(paste("NSE for validation period: ", round(NSE_validation, 3)))

## [1] "NSE for validation period:  0.858"

9.3 Compare the performance with the CemaNeige GR6J model

# Calculate NSE for GR6J model
OutputsCrit3 <- ErrorCrit_NSE(InputsCrit = InputsCrit3, OutputsModel = OutputsModel3)

## Crit. NSE[Q] = 0.7892

NSE_GR6J <- OutputsCrit3$CritValue 

# Calculate NSE for TUWmodel
NSE_TUW <- NSE(as.numeric(sim2$q)[366:konec], preob1$Qmm[366:konec])

# Compare the performance of the two models
if (NSE_GR6J > NSE_TUW) {
  cat("The GR6J model performs better with NSE =", NSE_GR6J,
      "compared to TUWmodel with NSE =", NSE_TUW, "\n")
} else if (NSE_GR6J < NSE_TUW) {
  cat("The TUWmodel performs better with NSE =", NSE_TUW,
      "compared to GR6J with NSE =", NSE_GR6J, "\n")
} else {
  cat("Both models perform equally well with NSE =", NSE_GR6J, "\n")
}

## The TUWmodel performs better with NSE = 0.8579348 compared to GR6J with NSE = 0.7891693

The comparison of model performance based on Nash-Sutcliffe Efficiency (NSE) reveals that the TUWmodel outperforms the GR6J model. The TUWmodel achieved an NSE value of 0.858, compared to the GR6J model’s NSE value of 0.789. The superior performance of the TUWmodel can be attributed to its more sophisticated representation of hydrological processes, potentially including enhanced snowmelt modeling and soil moisture dynamics. This indicates that the TUWmodel provides a more precise simulation of observed streamflow, particularly in regions or conditions where such processes are significant.

10. STOCHASTIC CLIMATE SIMULATOR

In this section, we explore the use of stochastic simulation techniques to model precipitation. First, a stochastic precipitation simulator will be fitted to measured data using the GWEX package. Ten realizations of 13-year simulations will then be generated and compared to the measured data to evaluate the simulator’s performance. Finally, the simulated precipitation data will serve as input for the previously calibrated hydrological model (CemaNeige GR6J) to assess its performance under stochastic conditions.

10.1 Fit a Stochastic Precipitation and Air Temperature Simulator

# Transform Precipitation Data to Matrix
precip_matrix <- as.matrix(main_data[, 8:12])

# Define Precipitation Observations for GWEX
observed_precip <- GwexObs(variable = 'Prec', date = as.Date(main_data[, 1], format = "%Y.%m.%d"), 
                           obs = precip_matrix)

# Estimate Parameters for Precipitation Simulator
precip_params <- fitGwexModel(observed_precip, listOption = list(th = 0.5, nChainFit = 1000))

## [1] "Fit generator"

precip_params

## *** Class GWex *** 
## * type of variable = Prec
## * Version = v1.1.0
## 
## #### Options for the occurrence process ####
## threshold: 0.5 
## order of the Markov chain: 2 
## 
## #### Options for the amount process ####
## spatial dependence: Gaussian copula 
## apply a MAR(1) process: FALSE 
## apply a master Markov chain:

10.2 Simulate Realizations and Compare with Measured Data

# Simulate Precipitation (10 Realizations for 13 Years)
simulated_precip <- simGwexModel(precip_params, nb.rep = 10, d.start = as.Date("01012025", "%d%m%Y"), 
                                 d.end = as.Date("31122037", "%d%m%Y"))

## [1] "Generate scenarios"

# Summary of Simulated Precipitation
precip_summary <- summary(simulated_precip@sim[, 1, 1])
precip_summary

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   4.089   2.486 183.289

P1_Luce_summary <- summary(main_data$P1_Luce)
P1_Luce_summary

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   4.236   2.200 116.100

# Plot Simulated Precipitation for Station 1, Realization 1
plot(simulated_precip@date, simulated_precip@sim[,1,1], type = "l", col = "green", 
     main = "Simulated Precipitation Station 1 (P1_Luce)", xlab = "Date", 
     ylab = "Precipitation (mm)")

# Add Realization 2
lines(simulated_precip@date, simulated_precip@sim[,1,2], col = "purple")

legend("topright", legend = c("Realization 1", "Realization 2"), col = c("green", "purple"), lty = 1)

# P1_Luce measured yearly precipitation
P1_Luce_yearly <- main_data %>% 
  group_by(Year) %>% 
  summarise(P1_Luce_measured = sum(P1_Luce))
P1_Luce_yearly

## # A tibble: 13 × 2
##     Year P1_Luce_measured
##    <dbl>            <dbl>
##  1  2010            1726.
##  2  2011            1338.
##  3  2012            1712.
##  4  2013            1538.
##  5  2014            2088.
##  6  2015            1261 
##  7  2016            1610.
##  8  2017            1652.
##  9  2018            1436.
## 10  2019            1553 
## 11  2020            1503.
## 12  2021            1554.
## 13  2022            1139.

# Dataframe of dates and years
simulated_yearly <- data.frame(Date = simulated_precip@date, 
                               Year = format(simulated_precip@date, "%Y"))

# Add yearly precipitation for each realization
for (i in 1:10) { 
  simulated_yearly[[paste0("Realization_", i)]] <- simulated_precip@sim[, 1, i] 
  }

# Calculate yearly totals for each realization
simulated_combined <- simulated_yearly %>% 
  group_by(Year) %>% 
  summarise(across(starts_with("Realization"), ~ sum(.x, na.rm = TRUE)))

simulated_combined

## # A tibble: 13 × 11
##    Year  Realization_1 Realization_2 Realization_3 Realization_4 Realization_5
##    <chr>         <dbl>         <dbl>         <dbl>         <dbl>         <dbl>
##  1 2025          1207.         1623.         1736.         1846.         1483.
##  2 2026          1541.         1417.         1396.         1534.         1476.
##  3 2027          1810.         1613.         1650.         1752.         1652.
##  4 2028          1431.         1558.         1804.         1185.         1772.
##  5 2029          1256.         1774.         1728.         1343.         1429.
##  6 2030          1401.         1410.         1427.         1524.         1562.
##  7 2031          1495.         1792.         1685.         1455.         1734.
##  8 2032          1590.         1108.         1183.         1577.         1627.
##  9 2033          1226.         1349.         1787.         2001.         1268.
## 10 2034          1658.         1241.         1922.         1388.         1082.
## 11 2035          1275.         1258.         1308.         1480.         1672.
## 12 2036          1964.         1550.         1705.         1535.         1620.
## 13 2037          1562.         1613.         1198.         1757.         1291.
## # ℹ 5 more variables: Realization_6 <dbl>, Realization_7 <dbl>,
## #   Realization_8 <dbl>, Realization_9 <dbl>, Realization_10 <dbl>

# Convert Year to numeric
P1_Luce_yearly <- P1_Luce_yearly %>% mutate(Year = as.numeric(Year))
simulated_combined <- simulated_combined %>% mutate(Year = as.numeric(Year))

# Set up the plot for measured precipitation
plot(P1_Luce_yearly$Year, P1_Luce_yearly$P1_Luce_measured, type = "o", col = "blue", pch = 16, lwd = 2, 
     xlab = "Year", ylab = "Total Precipitation (mm)", main = "Measured vs Simulated Yearly Precipitation(P1_Luce)", 
     ylim = c(500, 2500))

# Add simulated precipitation for each realization
for (i in 1:10) { 
  lines(P1_Luce_yearly$Year, simulated_combined[[paste0("Realization_", i)]],
  col = "red", lty = 2) 
}

# Add a legend
legend("topright", legend = c("Measured", "Simulated (Realizations)"), 
       col = c("blue", "red"), lty = c(1, 2), pch = c(16, NA), lwd = c(2, 1))

10.3 Use Simulated Data for Hydrological Models

for (i in 1:10) {
  calibration_end1 <- which(format(as.Date(simulated_precip@date), format = "%Y-%m-%d") == "2031-12-31")
  Ind_Run2 <- seq(366, calibration_end1)
  
  InputsModel3 <- CreateInputsModel(FUN_MOD = RunModel_CemaNeigeGR6J, 
                                    DatesR = as.POSIXct(simulated_precip@date, format = "%Y-%m-%d"), 
                                    Precip = simulated_precip@sim[, 1, i], PotEvap = model_data[, 5], 
                                    TempMean = model_data[, 6], ZInputs = median(hypso), 
                                    HypsoData = hypso, NLayers = 5)
  
  RunOptions3 <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR6J, InputsModel = InputsModel3, 
                                  IndPeriod_Run = Ind_Run2, IndPeriod_WarmUp = 1:365, 
                                  MeanAnSolidPrecip = c(70, 130, 180, 250, 350))
  
  OutputsModel3 <- RunModel_CemaNeigeGR6J(InputsModel = InputsModel3, RunOptions = RunOptions3, 
                                          Param = ParamGR6JCemaNeige)
  
  plot(OutputsModel3)
}

## input series were successfully created on 5 elevation layers for use by CemaNeige

## Warning in plot.OutputsModel(OutputsModel3): the "Error" and "CorQQ" plot(s)
## cannot be drawn if there is no 'Qobs'

## input series were successfully created on 5 elevation layers for use by CemaNeige

## Warning in plot.OutputsModel(OutputsModel3): the "Error" and "CorQQ" plot(s)
## cannot be drawn if there is no 'Qobs'

## input series were successfully created on 5 elevation layers for use by CemaNeige

## Warning in plot.OutputsModel(OutputsModel3): the "Error" and "CorQQ" plot(s)
## cannot be drawn if there is no 'Qobs'

## input series were successfully created on 5 elevation layers for use by CemaNeige

## Warning in plot.OutputsModel(OutputsModel3): the "Error" and "CorQQ" plot(s)
## cannot be drawn if there is no 'Qobs'

## input series were successfully created on 5 elevation layers for use by CemaNeige

## Warning in plot.OutputsModel(OutputsModel3): the "Error" and "CorQQ" plot(s)
## cannot be drawn if there is no 'Qobs'

## input series were successfully created on 5 elevation layers for use by CemaNeige

## Warning in plot.OutputsModel(OutputsModel3): the "Error" and "CorQQ" plot(s)
## cannot be drawn if there is no 'Qobs'

## input series were successfully created on 5 elevation layers for use by CemaNeige

## Warning in plot.OutputsModel(OutputsModel3): the "Error" and "CorQQ" plot(s)
## cannot be drawn if there is no 'Qobs'

## input series were successfully created on 5 elevation layers for use by CemaNeige

## Warning in plot.OutputsModel(OutputsModel3): the "Error" and "CorQQ" plot(s)
## cannot be drawn if there is no 'Qobs'

## input series were successfully created on 5 elevation layers for use by CemaNeige

## Warning in plot.OutputsModel(OutputsModel3): the "Error" and "CorQQ" plot(s)
## cannot be drawn if there is no 'Qobs'

## input series were successfully created on 5 elevation layers for use by CemaNeige

## Warning in plot.OutputsModel(OutputsModel3): the "Error" and "CorQQ" plot(s)
## cannot be drawn if there is no 'Qobs'

12. MODIS AND ERA5 DATA ANALYSIS

In this section, we analyze the surface air temperature data from MODIS for the Savinja_VSirje catchment over a period of at least two years. The temperature characteristics are examined and compared with measured discharge data to identify potential correlations and impacts on hydrological behavior. Additionally, ERA5 data for the entire available period is downloaded and analyzed. A comparison between ERA5 spatial precipitation data and gauge-based areal precipitation measurements is conducted, with plots provided to visually compare the datasets and highlight any differences or similarities.

12.1 Download and Process MODIS Surface Temperature Data

# List all available MODIS products
head(mt_products())

##        product
## 1       Daymet
## 2 ECO4ESIPTJPL
## 3      ECO4WUE
## 4       GEDI03
## 5     GEDI04_B
## 6      MCD12Q1
##                                                                          description
## 1 Daily Surface Weather Data (Daymet) on a 1-km Grid for North America, Version 4 R1
## 2               ECOSTRESS Evaporative Stress Index PT-JPL (ESI) Daily L4 Global 70 m
## 3                          ECOSTRESS Water Use Efficiency (WUE) Daily L4 Global 70 m
## 4                GEDI Gridded Land Surface Metrics (LSM) L3 1km EASE-Grid, Version 2
## 5     GEDI Gridded Aboveground Biomass Density (AGBD) L4B 1km EASE-Grid, Version 2.1
## 6              MODIS/Terra+Aqua Land Cover Type (LC) Yearly L3 Global 500 m SIN Grid
##   frequency resolution_meters
## 1     1 day              1000
## 2    Varies                70
## 3    Varies                70
## 4  One time              1000
## 5  One time              1000
## 6    1 year               500

# Check available bands for MODIS product MYD11A2
head(mt_bands("MYD11A2"))

##               band                          description valid_range fill_value
## 1   Clear_sky_days               Day clear-sky coverage    1 to 255          0
## 2 Clear_sky_nights             Night clear-sky coverage    1 to 255          0
## 3    Day_view_angl View zenith angle of day observation    0 to 130        255
## 4    Day_view_time        Local time of day observation    0 to 240        255
## 5          Emis_31                   Band 31 emissivity    1 to 255          0
## 6          Emis_32                   Band 32 emissivity    1 to 255          0
##    units scale_factor add_offset
## 1   <NA>         <NA>       <NA>
## 2   <NA>         <NA>       <NA>
## 3 degree            1        -65
## 4    hrs          0.1          0
## 5   <NA>        0.002       0.49
## 6   <NA>        0.002       0.49

# Download MODIS surface temperature data
modis_temp <- mt_subset(
  product = "MYD11A2",
  lat = 46.0959,
  lon = 15.1874,
  band = "LST_Day_1km",
  start = "2021-01-01",
  end = "2022-12-31",
  km_lr = 80,
  km_ab = 80,
  site_name = "Savinja-VelikoSirje",
  internal = TRUE,
  progress = FALSE
)

# Convert data to terra format
modis_temp_raster <- mt_to_terra(df = modis_temp)

# Plot surface temperature data (Kelvin)
plot(modis_temp_raster)

# Project the first layer (Jan 1, 2021)
selected_layer <- project(modis_temp_raster[[1]], crspost)
plot(selected_layer)
lines(catchment_area, col = "red")

# Extract temperature data for each date and convert to Celsius
temperature_data <- data.frame()

for (i in 1:nlyr(modis_temp_raster)) {
  projected_layer <- project(modis_temp_raster[[i]], crspost)
  projected_layer <- subst(projected_layer, 0, NA)
  temp_kelvin <- extract(projected_layer, catchment_area, mean, rm)[2]
  
  if (!is.na(temp_kelvin)) {
    temp_celsius <- temp_kelvin - 273.15
    new_row <- data.frame(Date = as.Date(names(modis_temp_raster)[i]), Temp_Celsius = temp_celsius, Temp_Kelvin = temp_kelvin)
    colnames(new_row) <- c("Date", "Temp_Celsius", "Temp_Kelvin")
    temperature_data <- rbind(temperature_data, new_row)
  }
}

head(temperature_data)

##         Date Temp_Celsius Temp_Kelvin
## 1 2021-02-10    0.7508639    273.9009
## 2 2021-02-18   13.5209944    286.6710
## 3 2021-02-26   13.3185523    286.4686
## 4 2021-03-06    7.8965254    281.0465
## 5 2021-03-22   16.3787452    289.5287
## 6 2021-03-30   21.1023461    294.2523

12.2 Comparison of MODIS and Measured Temperature Data

# Filter measured air temperature data
measured_temp <- main_data[main_data$Date >= as.Date("2021-01-01"), c("Date", "Airtemp")]

# Plot comparison between measured and MODIS temperature data
plot(
  as.Date(measured_temp$Date), measured_temp$Airtemp, 
  type = "l", col = "blue", xlab = "Date", ylab = "Temperature (°C)", 
  main = "Temperature Comparison: Measured vs. MODIS",
  ylim = c(-10, 35)
  
)
lines(as.Date(temperature_data$Date), temperature_data$Temp_Celsius, col = "red", type = "o", pch = 16)

legend(
  "topright", 
  legend = c("Measured Data", "MODIS Data"), 
  col = c("blue", "red"), 
  lty = 1, 
  pch = c(NA, 16)
)

12.3 Download and Process ERA5 Precipitation Data

# Load ERA5 precipitation data (yearly, land domain)
# download_data("era5", getwd(), timestep = "yearly", domain = "land")

era5_precip <- rast("era5_tp_mm_land_195901_202112_025_yearly.nc")

# View basic properties and plot the first layer
era5_precip

## class       : SpatRaster 
## dimensions  : 720, 1440, 63  (nrow, ncol, nlyr)
## resolution  : 0.25, 0.25  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
## source      : era5_tp_mm_land_195901_202112_025_yearly.nc 
## varname     : tp (Total precipitation) 
## names       : tp_1, tp_2, tp_3, tp_4, tp_5, tp_6, ... 
## unit        :   mm,   mm,   mm,   mm,   mm,   mm, ... 
## time        : 1959-01-01 to 2021-01-01 UTC

plot(era5_precip[[1]])

# Project catchment area to match ERA5 data
catchment_projected <- project(catchment_area, crs(era5_precip))

# Plot ERA5 data with catchment boundary overlay
plot(era5_precip[["tp_1"]], ext = c(14, 16, 45.5, 47))
lines(catchment_projected, col = "red")

# Extract grid cells within the catchment area
extracted_cells <- extract(era5_precip[[1]], catchment_projected, xy = TRUE, touches = TRUE, weights = TRUE)
head(extracted_cells)

##   ID     tp_1      x      y  weight
## 1  1 1236.111 14.625 46.375 0.35920
## 2  1 1234.035 14.875 46.375 0.71614
## 3  1 1278.850 15.125 46.375 0.75436
## 4  1 1233.905 15.375 46.375 0.43063
## 5  1 1129.915 15.625 46.375 0.00890
## 6  1 1233.344 14.875 46.125 0.16533

# Calculate weighted precipitation
weighted_precip <- sum(extracted_cells$tp_1 * extracted_cells$weight, na.rm = TRUE)
weighted_precip

## [1] 4197.993

annual_precip <- rep(NA, dim(era5_precip)[3])

# Compute average precipitation for each year
for (i in 1:dim(era5_precip)[3]) {
  values <- extract(era5_precip[[i]], catchment_projected)[[2]]
  annual_precip[i] <- mean(values, na.rm = TRUE)
}

# Plot annual ERA5 precipitation
time_series <- as.Date(time(era5_precip))
plot(
  time_series, annual_precip, type = "o", 
  xlab = "Year", ylab = "Annual Precipitation (mm)", 
  main = "Annual ERA5 Precipitation for Savinja Catchment"
)

12.4 Comparison of ERA5 and Areal Precipitation Data

# Aggregate daily Areal Precipitation to annual totals
areal_precip_annual <- main_data %>%
  group_by(Year) %>%
  summarise(Areal_Precip = sum(Areal_Precip, na.rm = TRUE))

# Convert Year to Date format for plotting
areal_precip_annual$Year <- as.Date(paste0(areal_precip_annual$Year, "-01-01"))

# Filter ERA5 data from 2010 onward
time_vector <- as.Date(time(era5_precip))
start_index <- which(time_vector == as.Date("2010-01-01"))
time_vector_filtered <- time_vector[start_index:length(time_vector)]
annual_precip_filtered <- annual_precip[start_index:length(annual_precip)]

# Plot Areal and ERA5 Precipitation Comparison
plot(
  areal_precip_annual$Year, areal_precip_annual$Areal_Precip, type = "o", col = "blue",
  xlab = "Year", ylab = "Precipitation (mm)", 
  main = "Comparison of Areal and ERA5 Precipitation",
  xlim = as.Date(c("2010-01-01", "2021-12-31")),
  ylim = c(900, 1700),
  xaxt = "n"
)

# Add ERA5 precipitation data
lines(time_vector_filtered, annual_precip_filtered, col = "red", lty = 1, type = "o")
axis(
  1,
  at = seq(as.Date("2010-01-01"), as.Date("2022-01-01"), by = "1 year"),
  labels = format(seq(as.Date("2010-01-01"), as.Date("2022-01-01"), by = "1 year"), "%Y")
)

# Add legend
legend(
  "topright", legend = c("Areal Precipitation", "ERA5 Precipitation"),
  col = c("blue", "red"), lty = 1, pch = c(16, 16)
)